def test_abu_evolution(self):
        from NuGridPy import ppn, utils
        import matplotlib
        matplotlib.use('agg')
        import matplotlib.pylab as mpy
        import os

        # Perform tests within temporary directory
        with TemporaryDirectory() as tdir:
            # wget the data for a ppn run from the CADC VOspace
            os.system("wget -q --content-disposition --directory '" + tdir +  "' "\
                          + "'http://www.canfar.phys.uvic.ca/vospace/synctrans?TARGET="\
                          + "vos%3A%2F%2Fcadc.nrc.ca%21vospace%2Fnugrid%2Fdata%2Fprojects%2Fppn%2Fexamples%2F"\
                          + "ppn_Hburn_simple%2Fx-time.dat&DIRECTION=pullFromVoSpace&PROTOCOL"\
                          + "=ivo%3A%2F%2Fivoa.net%2Fvospace%2Fcore%23httpget'")

            #nugrid_dir= os.path.dirname(os.path.dirname(ppn.__file__))
            #NuPPN_dir= nugrid_dir + "/NuPPN"
            #test_data_dir= NuPPN_dir + "/examples/ppn_Hburn_simple/RUN_MASTER"

            symbs=utils.symbol_list('lines2')
            x=ppn.xtime(tdir)
            specs=['PROT','HE  4','C  12','N  14','O  16']
            i=0
            for spec in specs:
                x.plot('time',spec,logy=True,logx=True,shape=utils.linestyle(i)[0],show=False,title='')
                i += 1
            mpy.ylim(-5,0.2)
            mpy.legend(loc=0)
            mpy.xlabel('$\log t / \mathrm{min}$')
            mpy.ylabel('$\log X \mathrm{[mass fraction]}$')
            abu_evol_file = 'abu_evolution.png'
            mpy.savefig(abu_evol_file)
            self.assertTrue(os.path.exists(abu_evol_file))

def plot_svc(X, y, mysvc, bounds=None, grid=50):
    if bounds is None:
        xmin = np.min(X[:, 0], 0)
        xmax = np.max(X[:, 0], 0)
        ymin = np.min(X[:, 1], 0)
        ymax = np.max(X[:, 1], 0)
    else:
        xmin, ymin = bounds[0], bounds[0]
        xmax, ymax = bounds[1], bounds[1]
    aspect_ratio = (xmax - xmin) / (ymax - ymin)
    xgrid, ygrid = np.meshgrid(np.linspace(xmin, xmax, grid),
                              np.linspace(ymin, ymax, grid))
    plt.gca(aspect=aspect_ratio)
    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.xticks([])
    plt.yticks([])
    plt.hold(True)
    plt.plot(X[y == 1, 0], X[y == 1, 1], 'bo')
    plt.plot(X[y == -1, 0], X[y == -1, 1], 'ro')
    
    box_xy = np.append(xgrid.reshape(xgrid.size, 1), ygrid.reshape(ygrid.size, 1), 1)
    if mysvc is not None:
        scores = mysvc.decision_function(box_xy)
    else:
        print 'You must have a valid SVC object.'
        return None;
    
    CS=plt.contourf(xgrid, ygrid, scores.reshape(xgrid.shape), alpha=0.5, cmap='jet_r')
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[0], colors='k', linestyles='solid', linewidths=1.5)
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[-1,1], colors='k', linestyles='dashed', linewidths=1)
    plt.plot(mysvc.support_vectors_[:,0], mysvc.support_vectors_[:,1], 'ko', markerfacecolor='none', markersize=10)
    CB = plt.colorbar(CS)

def study_multiband_planck(quick=True):
    savename = datadir+'cl_multiband.pkl'
    bands = [100, 143, 217, 'mb']
    if quick: cl = pickle.load(open(savename,'r'))
    else:
        cl = {}
        mask = load_planck_mask()
        mask_factor = np.mean(mask**2.)
        for band in bands:
            this_map = load_planck_data(band)
            this_cl = hp.anafast(this_map*mask, lmax=lmax)/mask_factor
            cl[band] = this_cl
        pickle.dump(cl, open(savename,'w'))


    cl_theory = {}
    pl.clf()
    
    for band in bands:
        l_theory, cl_theory[band] = get_cl_theory(band)
        this_cl = cl[band]
        pl.plot(this_cl/cl_theory[band])
        
    pl.legend(bands)
    pl.plot([0,4000],[1,1],'k--')
    pl.ylim(.7,1.3)
    pl.ylabel('data/theory')

def plot_part2(filename):
	"""
	Plots the result of count ones test
	"""
	fig1 = pl.figure()
	iterations, runtimes, fvals = extract(filename)
	algos = ["SA", "GA", "MIMIC"]
	iters_sa, iters_ga, iters_mimic = [np.array(iterations[a]) for a in algos]
	runtime_sa, runtime_ga, runtime_mimic = [np.array(runtimes[a]) for a in algos]
	fvals_sa, fvals_ga, fvals_mimic = [np.array(fvals[a]) for a in algos]

	plotfunc = getattr(pl, "loglog")
	plotfunc(runtime_sa, fvals_sa, "bs", mew=0)
	plotfunc(runtime_ga, fvals_ga, "gs", mew=0)
	plotfunc(runtime_mimic, fvals_mimic, "rs", mew=0)

	# plotfunc(iters_sa, fvals_sa/(runtime_sa * iters_sa), "bs", mew=0)
	# plotfunc(iters_ga, fvals_ga/(runtime_ga * iters_ga), "gs", mew=0)
	# plotfunc(iters_mimic, fvals_mimic/(runtime_mimic * iters_mimic), "rs", mew=0)

	pl.xlabel("Runtime (seconds)")
	pl.ylabel("Objective function value")
	pl.ylim([min(fvals_sa) / 2, max(fvals_mimic) * 2])
	pl.legend(["SA", "GA", "MIMIC"], loc=4)

	pl.savefig(filename.replace(".csv", ".png"), bbox_inches="tight") 

    def visualization2(self, sp_to_vis=None):
        if sp_to_vis:
            species_ready = list(set(sp_to_vis).intersection(self.all_sp_signatures.keys()))
        else:
            raise Exception('list of driver species must be defined')

        if not species_ready:
            raise Exception('None of the input species is a driver')

        for sp in species_ready:
            # Setting up figure
            plt.figure()
            plt.subplot(313)

            mon_val = OrderedDict()
            signature = self.all_sp_signatures[sp]
            for idx, mon in enumerate(list(set(signature))):
                if mon[0] == 'C':
                    mon_val[self.all_comb[sp][mon] + (-1,)] = idx
                else:
                    mon_val[self.all_comb[sp][mon]] = idx

            mon_rep = [0] * len(signature)
            for i, m in enumerate(signature):
                if m[0] == 'C':
                    mon_rep[i] = mon_val[self.all_comb[sp][m] + (-1,)]
                else:
                    mon_rep[i] = mon_val[self.all_comb[sp][m]]
            # mon_rep = [mon_val[self.all_comb[sp][m]] for m in signature]

            y_pos = numpy.arange(len(mon_val.keys()))
            plt.scatter(self.tspan[1:], mon_rep)
            plt.yticks(y_pos, mon_val.keys())
            plt.ylabel('Monomials', fontsize=16)
            plt.xlabel('Time(s)', fontsize=16)
            plt.xlim(0, self.tspan[-1])
            plt.ylim(0, max(y_pos))

            plt.subplot(312)

            for name in self.model.odes[sp].as_coefficients_dict():
                mon = name
                mon = mon.subs(self.param_values)
                var_to_study = [atom for atom in mon.atoms(sympy.Symbol)]
                arg_f1 = [numpy.maximum(self.mach_eps, self.y[str(va)][1:]) for va in var_to_study]
                f1 = sympy.lambdify(var_to_study, mon)
                mon_values = f1(*arg_f1)
                mon_name = str(name).partition('__')[2]
                plt.plot(self.tspan[1:], mon_values, label=mon_name)
            plt.ylabel('Rate(m/sec)', fontsize=16)
            plt.legend(bbox_to_anchor=(-0.1, 0.85), loc='upper right', ncol=1)

            plt.subplot(311)
            plt.plot(self.tspan[1:], self.y['__s%d' % sp][1:], label=parse_name(self.model.species[sp]))
            plt.ylabel('Molecules', fontsize=16)
            plt.legend(bbox_to_anchor=(-0.15, 0.85), loc='upper right', ncol=1)
            plt.suptitle('Tropicalization' + ' ' + str(self.model.species[sp]))

            # plt.show()
            plt.savefig('s%d' % sp + '.png', bbox_inches='tight', dpi=400)

def plot_prob_effector(sens, fpr, xmax=1, baserate=0.1):
    """Plots a line graph of P(effector|positive test) against
    the baserate of effectors in the input set to the classifier.
        
    The baserate argument draws an annotation arrow
    indicating P(pos|+ve) at that baserate
    """
    assert 0.1 <= xmax <= 1, "Max x axis value must be in range [0,1]"
    assert 0.01 <= baserate <= 1, "Baserate annotation must be in range [0,1]"
    baserates = pylab.arange(0, 1.05, xmax * 0.005)  
    probs = [p_correct_given_pos(sens, fpr, b) for b in baserates]
    pylab.plot(baserates, probs, 'r')
    pylab.title("P(eff|pos) vs baserate; sens: %.2f, fpr: %.2f" % (sens, fpr))
    pylab.ylabel("P(effector|positive)")
    pylab.xlabel("effector baserate")
    pylab.xlim(0, xmax)
    pylab.ylim(0, 1)
    # Add annotation arrow
    xpos, ypos = (baserate, p_correct_given_pos(sens, fpr, baserate))
    if baserate < xmax:
        if xpos > 0.7 * xmax:
            xtextpos = 0.05 * xmax
        else:
            xtextpos = xpos + (xmax-xpos)/5.
        if ypos > 0.5:
            ytextpos = ypos - 0.05
        else:
            ytextpos = ypos + 0.05
        pylab.annotate('baserate: %.2f, P(pos|+ve): %.3f' % (xpos, ypos), 
                       xy=(xpos, ypos), 
                       xytext=(xtextpos, ytextpos),
                       arrowprops=dict(facecolor='black', shrink=0.05))
    else:
        pylab.text(0.05 * xmax, 0.95, 'baserate: %.2f, P(pos|+ve): %.3f' %
                   (xpos, ypos))

def plotFirstTacROC(dataset):
    import matplotlib.pylab as plt
    from os.path import join
    from src.utils import PROJECT_DIR

    plt.figure(figsize=(6, 6))
    time_sampler = TimeSerieSampler(n_time_points=12)
    evaluator = Evaluator()
    time_series_idx = 0
    methods = {
        "cross_correlation": "Cross corr.   ",
        "kendall": "Kendall        ",
        "symbol_mutual": "Symbol MI    ",
        "symbol_similarity": "Symbol sim.",
    }
    for method in methods:
        print method
        predictor = SingleSeriesPredictor(good_methods[method], time_sampler)
        prediction = predictor.predictAllInstancesCombined(dataset, time_series_idx)
        roc_auc, fpr, tpr = evaluator.evaluate(prediction)
        plt.plot(fpr, tpr, label=methods[method] + " (auc = %0.3f)" % roc_auc)
    plt.legend(loc="lower right")
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel("False Positive Rate")
    plt.ylabel("True Positive Rate")
    plt.grid()
    plt.savefig(join(PROJECT_DIR, "output", "firstTACROC.pdf"))

    def plotClusters(self, clusters, Size=(11, 8)):
        """ Plots events belong to five different clusters.

        **Parameters**

        clusters : int (array or list)
            The index of the cluster from which the events will be plotted
        """
        fig = plt.figure(figsize=Size)
        fig.subplots_adjust(wspace=.3, hspace=.3)

        ax = fig.add_subplot(511)
        self.plot_event(self.evts[self.goodEvts, :]
                        [np.array(clusters) == 0, :])
        plt.ylim([-15, 20])

        ax = fig.add_subplot(512)
        self.plot_event(self.evts[self.goodEvts, :]
                        [np.array(clusters) == 1, :])
        ax.set_ylim([-15, 20])

        ax = fig.add_subplot(513)
        self.plot_event(self.evts[self.goodEvts, :]
                        [np.array(clusters) == 2, :])
        ax.set_ylim([-15, 20])

        ax = fig.add_subplot(514)
        self.plot_event(self.evts[self.goodEvts, :]
                        [np.array(clusters) == 3, :])
        ax.set_ylim([-15, 20])

        ax = fig.add_subplot(515)
        self.plot_event(self.evts[self.goodEvts, :]
                        [np.array(clusters) == 4, :])
        ax.set_ylim([-15, 20])

def km_emp_mean(df_pca,krange,empcol,crtcol):
    '''
    apply kmeans to pca-processed dataframe, with a range of k. Then check with k yields clearest employed rate
    and correct rate.
    plot total number of k against good cluster ratio
    :param df_pca: numpy array
    :param krange: int, range of k
    :param empcol: panda series
    :param crtcol: panda series
    :raise: string, data frame, and plot
    '''
    lis = []
    df = pd.concat([empcol,crtcol],axis=1)
    for k in xrange(4,krange):
        km = KMeans(n_clusters=k,random_state=0)
        km.fit(df_pca)
        df['cluster'] = km.labels_
        temp = df.groupby('cluster').agg(np.mean)
        res = temp[(temp['employed'] >0.7) | (temp['employed'] < 0.4) &(temp['correct'] >0.7) ]
        print "{} clusters".format(k)
        print res
        print "---"
        print "{} out of {} clusters split the target ideally. Good cluster rate: {}".format(res.shape[0],k,res.shape[0]/k)
        print '*'*20
        lis.append(res.shape[0]/k)
    plt.plot(range(4,krange),lis,lw = 2)
    plt.xlabel('K')
    plt.ylabel('good cluster rate')
    plt.ylim(0.5,1.2)

def plot_eye(Nodes,axes = None):
    """
    
    Create a movie of eye growth. To be used with EyeGrowthFunc

    :param Nodes: structure containing nodes
    :type Nodes: struct
    :param INTSTEP: time step used for integration
    :type INTSTEP: int

    :returns: plot handle for Node plot.  Used to update during for loop.

    .. note::
       Called in EyeGrowthFunc
    """
    
    
    #set plotting parameters:
    if axes == None:
        fig = plt.figure(figsize=(10, 8))
        axes = fig.add_subplot(111,aspect='equal')
        plt.xlim([-13, 13])
        plt.ylim([-13, 13])
        
    axes.plot(np.r_[ Nodes['x'][0],Nodes['x'][0,0] ] * Nodes['radius'], 
            np.r_[ Nodes['y'][0], Nodes['y'][0,0] ] * Nodes['radius'], 
            '-ok', markerfacecolor = 'k',linewidth = 4, markersize = 10)
                
    axes = pf.TufteAxis(axes,['left','bottom'])
    #axes.set_axis_bgcolor('w')


    return axes

def plotParamPath(chain, OBJi,
                  save_opt=None):
    import matplotlib.pylab as plt
    import variable_house as vHouse
    nWalkers = np.shape(chain)[0]
    nSteps = np.shape(chain)[1]
    nDim = np.shape(chain)[2]
    stepidx = np.arange(1,nSteps+1)
    
    for param_i in range(0,nDim-vHouse.nDram):
        for walker_i in range(0,nWalkers):
                plt.plot(stepidx, chain[walker_i,:,param_i],
                         color='k',
                         alpha= 0.1)
        
        plt.xlabel('STEP INDEX')
        plt.ylabel(vHouse.nameL[param_i])
        plt.title('param value of walkers after trimming')
        #show entire range of possible parameter values
        plt.ylim([np.min(vHouse.paramRanges[param_i,:]),
                  np.max(vHouse.paramRanges[param_i,:])]) 
        if save_opt != None:
            plt.savefig(save_opt+'param'+str(param_i)+'path.png')
            plt.close('all')
        else:
            plt.show()   

def demo():
    '''
    Load and plot a few CIB spectra.
    '''

    # define ell array.
    l = np.arange(100,4000)

    # get dictionary of CIBxCIB spectra.
    cl_cibcib = get_cl_cibcib(l)

    # plot
    import matplotlib.pylab as pl
    pl.ion()
    lw=2
    fs=18
    leg = []
    pl.clf()
    for band in ['857','545','353']:
        pl.semilogy(l, cl_cibcib['545',band],linewidth=lw)
        leg.append('545 x '+band)
    pl.xlabel(r'$\ell$',fontsize=fs)
    pl.ylabel(r'$C_\ell^{TT, CIB} [\mu K^2]$',fontsize=fs)
    pl.ylim(5e-2,6e3)
    pl.legend(leg, fontsize=fs)

def read_table(args):
    tblfilename = "bf_optimize_mavlink.h5"
    h5file = tb.open_file(tblfilename, mode = "r")
    # print h5file
    # a = h5file.root
    # print a
    # a = h5file.get_node(where = "/20150507-run1/params/_20150507155408")
    # print a
    # table = h5file.root.v1.evaluations
    table = h5file.root.v2.evaluations
    print "table", table
    # mse = [x["mse"] for x in table.iterrows() if x["alt_p"] < 20.]
    # mse = [x["mse"] for x in table.iterrows()]
    logdata = [x["timeseries"] for x in table.iterrows() if x["mse"] < 2000]
    alt_pid = [(x["alt_p"], x["alt_i"], x["alt_d"], x["vel_p"], x["vel_i"], x["vel_d"]) for x in table.iterrows() if x["mse"] < 1000]
    # alt_pid = [(x["alt_p"], x["alt_i"], x["alt_d"]) for x in table.iterrows() if x["alt_p"] == 17 and x["alt_i"] == 0.]
    print "alt_pid", alt_pid
    # print mse
    # pl.plot(mse)
    print len(logdata)
    for i in range(len(logdata)):
        pl.subplot(len(logdata), 1, i+1)
        pl.plot(logdata[i][:,1:3])
        pl.ylim((-300, 1000))
    pl.show()

def plot_lift_data(lift_data, with_ellipses=True):
    np.random.seed(42113)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    alpha = [l['fit']['alpha'] for l in lift_data.values()]
    alpha_error = [l['fit']['alpha_error'] for l in lift_data.values()]
    beta = [l['fit']['beta'] for l in lift_data.values()]
    beta_error = [l['fit']['beta_error'] for l in lift_data.values()]
    message_class = lift_data.keys()

    num = len(beta)
    beta_jitter = np.random.randn(num)
    np.random.seed(None)
    beta = np.array(beta) + beta_jitter*0.0

    ax.plot(beta, alpha, color='red', linestyle='', marker='o', markersize=10)
    if not with_ellipses:
        ax.errorbar(beta, alpha, xerr=beta_error, yerr=alpha_error, linestyle='')
    else:
        for x, y, xerr, yerr, in zip(beta, alpha, beta_error, alpha_error):
            width = 2*xerr
            height = 2*yerr
            ellipse = patches.Ellipse((x, y), width, height,
                                      angle=0.0, linewidth=2,
                                      fill=True, alpha=0.15, color='gray')
            ax.add_patch(ellipse)

    for a, b, c in zip(alpha, beta, message_class):
        ax.annotate(c, xy=(b, a), xytext=(b+2, a+.01), fontsize=17)
    plt.xlim(0, max(beta)+30)
    plt.ylim(0, 0.9)
    plt.xlabel('Duration (days)')
    plt.ylabel('Initial Lift')
    plt.show()

def plot_average(filenames, save_plot=True, show_plot=False, dpi=100):

    ''' Plot Signal average from a list of averaged files. '''

    fname = get_files_from_list(filenames)

    # plot averages
    pl.ioff()  # switch off (interactive) plot visualisation
    factor = 1e15
    for fnavg in fname:
        name = fnavg[0:len(fnavg) - 4]
        basename = os.path.splitext(os.path.basename(name))[0]
        print fnavg
        # mne.read_evokeds provides a list or a single evoked based on condition.
        # here we assume only one evoked is returned (requires further handling)
        avg = mne.read_evokeds(fnavg)[0]
        ymin, ymax = avg.data.min(), avg.data.max()
        ymin *= factor * 1.1
        ymax *= factor * 1.1
        fig = pl.figure(basename, figsize=(10, 8), dpi=100)
        pl.clf()
        pl.ylim([ymin, ymax])
        pl.xlim([avg.times.min(), avg.times.max()])
        pl.plot(avg.times, avg.data.T * factor, color='black')
        pl.title(basename)

        # save figure
        fnfig = os.path.splitext(fnavg)[0] + '.png'
        pl.savefig(fnfig, dpi=dpi)

    pl.ion()  # switch on (interactive) plot visualisation

def plot_runs(runs):
    """ Plot population evolutions
    """
    ts = range(len(runs[0]))
    cmap = plt.get_cmap('viridis')
    for i, r in enumerate(runs):
        mean, var = zip(*r)
        bm, cm = zip(*mean)
        bv, cv = zip(*var)

        color = cmap(float(i)/len(runs))

        plt.errorbar(ts, bm, fmt='-', yerr=bv, c=color)
        plt.errorbar(ts, cm, fmt='--', yerr=cv, c=color)

    plt.title('population evolution overview')
    plt.xlabel('time')
    plt.ylabel('value')

    plt.ylim((0, 1))

    plt.plot(0, 0, '-', c='black', label='benefit value')
    plt.plot(0, 0, '--', c='black', label='cost value')
    plt.legend(loc='best')

    plt.savefig('result.pdf')
    plt.show()

        def _on_button_press(event):
            if event.button != 1 or not event.inaxes:
                return
            lon, lat = m(event.xdata, event.ydata, inverse=True)
            # Convert to colat to ease indexing.
            colat = rotations.lat2colat(lat)

            x_range = (self.setup["physical_boundaries_x"][1] -
                self.setup["physical_boundaries_x"][0])
            x_frac = (colat - self.setup["physical_boundaries_x"][0]) / x_range
            x_index = int(((self.setup["boundaries_x"][1] -
                self.setup["boundaries_x"][0]) * x_frac) +
                self.setup["boundaries_x"][0])
            y_range = (self.setup["physical_boundaries_y"][1] -
                self.setup["physical_boundaries_y"][0])
            y_frac = (lon - self.setup["physical_boundaries_y"][0]) / y_range
            y_index = int(((self.setup["boundaries_y"][1] -
                self.setup["boundaries_y"][0]) * y_frac) +
                self.setup["boundaries_y"][0])

            plt.figure(1, figsize=(3, 8))
            depths = available_depths
            values = data[x_index, y_index, :]
            plt.plot(values, depths)
            plt.grid()
            plt.ylim(depths[-1], depths[0])
            plt.show()
            plt.close()
            plt.figure(0)

 def plot(x,y,field,filename,c=200):
     plt.figure()
     # define grid.
     xi = np.linspace(min(x),max(x),100)
     yi = np.linspace(min(y),max(y),100)
     # grid the data.
     si_lin = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     si_cub = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     print np.min(field)
     print np.max(field)
     plt.subplot(211)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_lin,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_lin,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Lineaarinen interpolointi')
     #plt.tight_layout()
     plt.subplot(212)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_cub,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_cub,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Kuubinen interpolointi')
     plt.savefig(filename)

def draw_lineplot(x, y, title="title", xlab="x", ylab="y", odir="", xlim=None, ylim=None, outfmt='eps'):

  if len(x) == 0 or len(y) == 0:
    return;
  #fi

  plt.cla();
  plt.plot(x, y, marker='x');
  plt.xlabel(xlab);
  plt.ylabel(ylab);
  plt.title(title);

  if xlim == None:
    xmin = min(x);
    xmax = max(x);
    xlim = [xmin, xmax];
  #fi

  if ylim == None:
    ymin = min(y);
    ymax = max(y);
    ylim = [ymin, ymax];
  #fi

  plt.xlim(xlim);
  plt.ylim(ylim);

  plt.savefig('%s%s.%s' % (odir + ('/' if odir else ""), '_'.join(title.split(None)), outfmt), format=outfmt);

  return '%s.%s' % ('_'.join(title.split(None)), outfmt), title;

def plot_roc(test_cat, plot_data, savefig=False):
    # Plot ROC curve
    plt.clf()
    results = []

    # calcualte and sort labels by roc_auc
    for method, method_results in plot_data.items():
        fpr, tpr, roc_auc = compute_roc(test_cat, method_results, method)
        label = "[%s] area = %0.2f" % (method, roc_auc)
        res = {"label": label, "fpr": fpr, "tpr": tpr, "roc_auc": roc_auc}
        results.append(res)
    results = sorted(results, key=lambda k: k['roc_auc'], reverse=True)

    # plot according to roc_auc ranking
    for r in results:
        plt.plot(r["fpr"], r["tpr"], label=r["label"])

    plt.title('Receiver Operating Characteristic Curve (ROC)')
    plt.legend(loc="lower right")
    plt.plot([0, 1], [0, 1], 'k--')

    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.0])

    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')

    if savefig:
        plt.savefig("classifiers_roc.png")

    plt.show()