def check_models(self):
        temp = np.logspace(0, np.log10(600))
        num = len(self.available_models())

        fig, ax = plt.subplots(1)
        self.plotting_colours(num, fig, ax, repeats=2)

        for author in self.available_models():
            Nc, Nv = self.update(temp=temp, author=author)
            # print Nc.shape, Nv.shape, temp.shape
            ax.plot(temp, Nc, '--')
            ax.plot(temp, Nv, '.', label=author)

        ax.loglog()
        leg1 = ax.legend(loc=0, title='colour legend')

        Nc, = ax.plot(np.inf, np.inf, 'k--', label='Nc')
        Nv, = ax.plot(np.inf, np.inf, 'k.', label='Nv')

        plt.legend([Nc, Nv], ['Nc', 'Nv'], loc=4, title='Line legend')
        plt.gca().add_artist(leg1)

        ax.set_xlabel('Temperature (K)')
        ax.set_ylabel('Density of states (cm$^{-3}$)')
        plt.show()

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 make_corr1d_fig(dosave=False):
    corr = make_corr_both_hemi()
    lw=2; fs=16
    pl.figure(1)#, figsize=(8, 7))
    pl.clf()
    pl.xlim(4,300)
    pl.ylim(-400,+500)    
    lambda_titles = [r'$20 < \lambda < 30$',
                     r'$30 < \lambda < 40$',
                     r'$\lambda > 40$']
    colors = ['blue','green','red']
    for i in range(3):
        corr1d, rcen = corr_1d_from_2d(corr[i])
        ipdb.set_trace()
        pl.semilogx(rcen, corr1d*rcen**2, lw=lw, color=colors[i])
        #pl.semilogx(rcen, corr1d*rcen**2, 'o', lw=lw, color=colors[i])
    pl.xlabel(r'$s (Mpc)$',fontsize=fs)
    pl.ylabel(r'$s^2 \xi_0(s)$', fontsize=fs)    
    pl.legend(lambda_titles, 'lower left', fontsize=fs+3)
    pl.plot([.1,10000],[0,0],'k--')
    s_bao = 149.28
    pl.plot([s_bao, s_bao],[-9e9,+9e9],'k--')
    pl.text(s_bao*1.03, 420, 'BAO scale')
    pl.text(s_bao*1.03, 370, '%0.1f Mpc'%s_bao)
    if dosave: pl.savefig('xi1d_3bin.pdf')

def cdf(x,colsym="",lab="",lw=4):
    """ plot the cumulative density function

    Parameters
    ----------

    x : np.array()
    colsym : string
    lab : string
    lw : int
        linewidth

    Examples
    --------

    >>> import numpy as np

    """
    rcParams['legend.fontsize']=20
    rcParams['font.size']=20

    x  = np.sort(x)
    n  = len(x)
    x2 = np.repeat(x, 2)
    y2 = np.hstack([0.0, repeat(np.arange(1,n) / float(n), 2), 1.0])
    plt.plot(x2,y2,colsym,label=lab,linewidth=lw)
    plt.grid('on')
    plt.legend(loc=2)
    plt.xlabel('Ranging Error[m]')
    plt.ylabel('Cumulative Probability')

def plotMassFunction(im, pm, outbase, mmin=9, mmax=13, mstep=0.05):
    """
    Make a comparison plot between the input mass function and the 
    predicted projected correlation function
    """
    plt.clf()

    nmbins = ( mmax - mmin ) / mstep
    mbins = np.logspace( mmin, mmax, nmbins )
    mcen = ( mbins[:-1] + mbins[1:] ) /2
    
    plt.xscale( 'log', nonposx = 'clip' )
    plt.yscale( 'log', nonposy = 'clip' )
    
    ic, e, p = plt.hist( im, mbins, label='Original Halos', alpha=0.5, normed = True)
    pc, e, p = plt.hist( pm, mbins, label='Added Halos', alpha=0.5, normed = True)
    
    plt.legend()
    plt.xlabel( r'$M_{vir}$' )
    plt.ylabel( r'$\frac{dN}{dM}$' )
    #plt.tight_layout()
    plt.savefig( outbase+'_mfcn.png' )
    
    mdtype = np.dtype( [ ('mcen', float), ('imcounts', float), ('pmcounts', float) ] )
    mf = np.ndarray( len(mcen), dtype = mdtype )
    mf[ 'mcen' ] = mcen
    mf[ 'imcounts' ] = ic
    mf[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_mfcn.fit', mf )

def _fig_density(sweight, surweight, pval, nlm):
    """
    Plot the histogram of sweight across the image
    and the thresholds implied by the surrogate model (surweight)
    """
    import matplotlib.pylab as mp
    # compute some thresholds
    nlm = nlm.astype('d')
    srweight = np.sum(surweight,1)
    srw = np.sort(srweight)
    nitem = np.size(srweight)
    thf = srw[int((1-min(pval,1))*nitem)]
    mnlm = max(1,nlm.mean())
    imin = min(nitem-1,int((1.-pval/mnlm)*nitem))
    
    thcf = srw[imin]
    h,c = np.histogram(sweight,100)
    I = h.sum()*(c[1]-c[0])
    h = h/I
    h0,c0 = np.histogram(srweight,100)
    I0 = h0.sum()*(c0[1]-c0[0])
    h0 = h0/I0
    mp.figure(1)
    mp.plot(c,h)
    mp.plot(c0,h0)
    mp.legend(('true histogram','surrogate histogram'))
    mp.plot([thf,thf],[0,0.8*h0.max()])
    mp.text(thf,0.8*h0.max(),'p<0.2, uncorrected')
    mp.plot([thcf,thcf],[0,0.5*h0.max()])
    mp.text(thcf,0.5*h0.max(),'p<0.05, corrected')
    mp.savefig('/tmp/histo_density.eps')
    mp.show()

    def _plot_nullclines(self, resolution):
        """
        Plot nullclines.

        Arguments
            resolution
                Resolution of plot
        """
        x_mesh, y_mesh, ode_x, ode_y = self._get_ode_values(resolution)

        plt.contour(
            x_mesh, y_mesh, ode_x,
            levels=[0], linewidths=2, colors='black')
        plt.contour(
            x_mesh, y_mesh, ode_y,
            levels=[0], linewidths=2, colors='black',
            linestyles='dashed')

        lblx = mlines.Line2D(
            [], [],
            color='black',
            marker='', markersize=15,
            label=r'$\dot\varphi_0=0$')
        lbly = mlines.Line2D(
            [], [],
            color='black', linestyle='dashed',
            marker='', markersize=15,
            label=r'$\dot\varphi_1=0$')
        plt.legend(handles=[lblx, lbly], loc='best')

def find_params():

    FRAMES =  np.arange(30)*100

    frame_images = organizedata.get_frames(ddir("bukowski_04.W2"), FRAMES)
    print "DONE READING DATA"

    CLUST_EPS = np.linspace(0, 0.5, 10)
    MIN_SAMPLES = [2, 3, 4, 5]
    MIN_DISTS = [2, 3, 4, 5, 6]
    THOLD = 240

    fracs_2 = np.zeros((len(CLUST_EPS), len(MIN_SAMPLES), len(MIN_DISTS)))

    for cei, CLUST_EP in enumerate(CLUST_EPS):
        for msi, MIN_SAMPLE in enumerate(MIN_SAMPLES):
            for mdi, MIN_DIST in enumerate(MIN_DISTS):
                print cei, msi, mdi
                numclusters = np.zeros(len(FRAMES))
                for fi, im in enumerate(frame_images):
                    centers = frame_clust_points(im, THOLD, MIN_DIST, 
                                                 CLUST_EP, MIN_SAMPLE)
                    # cluster centers
                    numclusters[fi] = len(centers)
                fracs_2[cei, msi, mdi] = float(np.sum(numclusters == 2))/len(numclusters)
    pylab.figure(figsize=(12, 8))
    for mdi, MIN_DIST in enumerate(MIN_DISTS):
        pylab.subplot(len(MIN_DISTS), 1, mdi+1)

        for msi in range(len(MIN_SAMPLES)):
            pylab.plot(CLUST_EPS, fracs_2[:, msi, mdi], label='%d' % MIN_SAMPLES[msi])
        pylab.title("min_dist= %3.2f" % MIN_DIST)
    pylab.legend()
    pylab.savefig('test.png', dpi=300)

def test():
	## Load files
    s = load_spectrum('ring28yael')
    w = linspace(1510e-9,1600e-9,len(s))
    
	## Process
    mins = find_minima(s)
    w_p = 1510e-9 + array(mins) * 90.e-9/len(w)
    ww = 2 * pi * 3e8/w_p   
    
	## Plot
    pl.plot(w,s)
    pl.plot(w_p,s[mins],'o')
    pl.show()
    
    beta2 = -1./(112e-6*2*pi)*diff(diff(ww))/(diff(ww)[:-1]**3)
    p = polyfit(w_p[1:-1], beta2, 1)
    
    savetxt('ring28yael-p.txt', w_p)
    
    pl.subplot(211)
    pl.plot(w,s)
    pl.plot(w_p,s[mins],'o')
    
    pl.subplot(212)
    pl.plot(w_p[1:-1]*1e6, beta2)
    pl.plot(w_p[1:-1]*1e6, p[1]+ p[0]*w_p[1:-1], label="q=%.2e"%p[0])
    pl.legend()
        
    pl.show()

	def behavioral_analysis(self):
		"""some analysis of the behavioral data, such as mean percept duration, 
		dominance ratio etc"""
		self.assert_data_intern()
		# only do anything if this is not a no report trial
		if 'RP' in self.file_alias:
			all_percepts_and_durations = [[],[]]
		else:
			all_percepts_and_durations = [[],[],[]]
		if not 'NR' in self.file_alias: #  and not 'RP' in self.file_alias
			for x in range(len(self.trial_indices)):
				if len(self.events) != 0:
					events_this_trial = self.events[(self.events['EL_timestamp'] > self.timestamps_pt[x][0]) & (self.events['EL_timestamp'] < self.timestamps_pt[x][-1])]
					for sc, scancode in enumerate(self.scancode_list):
						percept_start_indices = np.arange(len(events_this_trial))[np.array(events_this_trial['scancode'] == scancode)]
						percept_end_indices = percept_start_indices + 1
						
						# convert to times
						start_times = np.array(events_this_trial['EL_timestamp'])[percept_start_indices] - self.timestamps_pt[x,0]
						if len(start_times) > 0:
							if percept_end_indices[-1] == len(events_this_trial):
								end_times = np.array(events_this_trial['EL_timestamp'])[percept_end_indices[:-1]] - self.timestamps_pt[x,0]
								end_times = np.r_[end_times, len(self.from_zero_timepoints)]
							else:
								end_times = np.array(events_this_trial['EL_timestamp'])[percept_end_indices] - self.timestamps_pt[x,0]

							these_raw_event_times = np.array([start_times + self.timestamps_pt[x,0], end_times + self.timestamps_pt[x,0]]).T
							these_event_times = np.array([start_times, end_times]).T + x * self.trial_duration * self.sample_rate
							durations = np.diff(these_event_times, axis = -1)

							all_percepts_and_durations[sc].append(np.hstack((these_raw_event_times, these_event_times, durations)))

			self.all_percepts_and_durations = [np.vstack(apd) for apd in all_percepts_and_durations]

			# last element is duration, sum inclusive and exclusive of transitions
			total_percept_duration = np.concatenate([apd[:,-1] for apd in self.all_percepts_and_durations]).sum()
			total_percept_duration_excl = np.concatenate([apd[:,-1] for apd in [self.all_percepts_and_durations[0], self.all_percepts_and_durations[-1]]]).sum()

			self.ratio_transition = 1.0 - (total_percept_duration_excl / total_percept_duration)
			self.ratio_percept_red = self.all_percepts_and_durations[0][:,-1].sum() / total_percept_duration_excl

			self.red_durations = np.array([np.mean(self.all_percepts_and_durations[0][:,-1]), np.median(self.all_percepts_and_durations[0][:,-1])])
			self.green_durations = np.array([np.mean(self.all_percepts_and_durations[-1][:,-1]), np.median(self.all_percepts_and_durations[-1][:,-1])])
			self.transition_durations = np.array([np.mean(self.all_percepts_and_durations[1][:,-1]), np.median(self.all_percepts_and_durations[1][:,-1])])

			self.ratio_percept_red_durations = self.red_durations / (self.red_durations + self.green_durations)
			plot_mean_or_median = 0 # mean

			f = pl.figure(figsize = (8,4))
			s = f.add_subplot(111)
			for i in range(len(self.colors)):
				pl.hist(self.all_percepts_and_durations[i][:,-1], bins = 20, color = self.colors[i], histtype='step', lw = 3.0, alpha = 0.4, label = ['Red', 'Trans', 'Green'][i])
			pl.hist(np.concatenate([self.all_percepts_and_durations[0][:,-1], self.all_percepts_and_durations[-1][:,-1]]), bins = 20, color = 'k', histtype='step', lw = 3.0, alpha = 0.4, label = 'Percepts')
			pl.legend()
			s.set_xlabel('time [ms]')
			s.set_ylabel('count')
			sn.despine(offset=10)
			s.annotate("""ratio_transition: %1.2f, \nratio_percept_red: %1.2f, \nduration_red: %2.2f,\nduration_green: %2.2f, \nratio_percept_red_durations: %1.2f"""%(self.ratio_transition, self.ratio_percept_red, self.red_durations[plot_mean_or_median], self.green_durations[plot_mean_or_median], self.ratio_percept_red_durations[plot_mean_or_median]), (0.5,0.65), textcoords = 'figure fraction')
			pl.tight_layout()
			pl.savefig(os.path.join(self.analyzer.fig_dir, self.file_alias + '_dur_hist.pdf'))

 def sanity_PDMAna(self):
   import numpy
   import matplotlib.pylab as mpl
   from PyAstronomy.pyTiming import pyPDM
   
   # Create artificial data with frequency = 3,
   # period = 1/3
   x = numpy.arange(100) / 100.0
   y = numpy.sin(x*2.0*numpy.pi*3.0 + 1.7)
   
   # Get a ``scanner'', which defines the frequency interval to be checked.
   # Alternatively, also periods could be used instead of frequency.
   S = pyPDM.Scanner(minVal=0.5, maxVal=5.0, dVal=0.01, mode="frequency")
   
   # Carry out PDM analysis. Get frequency array
   # (f, note that it is frequency, because the scanner's
   # mode is ``frequency'') and associated Theta statistic (t).
   # Use 10 phase bins and 3 covers (= phase-shifted set of bins).
   P = pyPDM.PyPDM(x, y)
   f1, t1 = P.pdmEquiBinCover(10, 3, S)
   # For comparison, carry out PDM analysis using 10 bins equidistant
   # bins (no covers).
   f2, t2 = P.pdmEquiBin(10, S)
   
   
   # Show the result
   mpl.figure(facecolor='white')
   mpl.title("Result of PDM analysis")
   mpl.xlabel("Frequency")
   mpl.ylabel("Theta")
   mpl.plot(f1, t1, 'bp-')
   mpl.plot(f2, t2, 'gp-')
   mpl.legend(["pdmEquiBinCover", "pdmEquiBin"])

def plotRocCurves(file_legend):
	pylab.clf()
	pylab.figure(1)
	pylab.xlabel('1 - Specificity', fontsize=12)
	pylab.ylabel('Sensitivity', fontsize=12)
	pylab.title("Need for Referral")
	pylab.grid(True, which='both')
	pylab.xticks([i/10.0 for i in range(1,11)])
	pylab.yticks([i/10.0 for i in range(0,11)])
	pylab.tick_params(axis="both", labelsize=15)

	for file, legend in file_legend:
		points = open(file,"rb").readlines()
		x = [float(p.split()[0]) for p in points]
		y = [float(p.split()[1]) for p in points]
		dev = [float(p.split()[2]) for p in points]
		x = [0.0] + x
		y = [0.0] + y
		dev = [0.0] + dev
	
		auc = np.trapz(y, x) * 100
		aucDev = np.trapz(dev, x) * 100

		pylab.grid()
		pylab.errorbar(x, y, yerr = dev, fmt='-')
		pylab.plot(x, y, '-', linewidth = 1.5, label = legend + u" (AUC = {0:0.1f}% \xb1 {1:0.1f}%)".format(auc,aucDev))

	pylab.legend(loc = 4, borderaxespad=0.4, prop={'size':12})
	pylab.savefig("referral/referral-curves.pdf", format='pdf')

def is_stationary(ts, test_window):
    """
	This function checks whether the given TS is stationary. Can make it boolean, but lets just leave it
	for visualisation purposes. Not to be run once the numbers have been fixed.
	"""

    # Determine the rolling statistics (places like these compelled me to use Pandas and not numpy here)
    rol_mean = pd.rolling_mean(ts, window=test_window)
    rol_std = pd.rolling_std(ts, window=test_window)

    # Plot rolling statistics:
    orig = plt.plot(ts, color="blue", label="Original")
    mean = plt.plot(rol_mean, color="red", label="Rolling Mean")
    std = plt.plot(rol_std, color="black", label="Rolling Std")
    plt.legend(loc="best")
    plt.title("Rolling Mean & Standard Deviation")
    plt.show()

    # Perform the  Dickey-Fuller test: (Check documentation of fn for return params)
    print "Results of Dickey-Fuller Test:"
    dftest = adfuller(timeseries, autolag="AIC")
    dfoutput = pd.Series(dftest[0:4], index=["Test Statistic", "p-value", "#Lags Used", "Number of Observations Used"])
    for key, value in dftest[4].items():
        dfoutput["Critical Value (%s)" % key] = value
    print dfoutput

    def check_models(self):
        plt.figure('Bandgap narrowing')
        Na = np.logspace(12, 20)
        Nd = 0.
        dn = 1e14
        temp = 300.

        for author in self.available_models():
            BGN = self.update(Na=Na, Nd=Nd, nxc=dn,
                              author=author,
                              temp=temp)

            if not np.all(BGN == 0):
                plt.plot(Na, BGN, label=author)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'Bgn.csv')

        data = np.genfromtxt(test_file, delimiter=',', names=True)

        for name in data.dtype.names[1:]:
            plt.plot(
                data['N'], data[name], 'r--',
                label='PV-lighthouse\'s: ' + name)

        plt.semilogx()
        plt.xlabel('Doping (cm$^{-3}$)')
        plt.ylabel('Bandgap narrowing (K)')

        plt.legend(loc=0)

def simulationWithoutDrug(numViruses, maxPop, maxBirthProb, clearProb,
                          numTrials):
    """
    Run the simulation and plot the graph for problem 3 (no drugs are used,
    viruses do not have any drug resistance).    
    For each of numTrials trial, instantiates a patient, runs a simulation
    for 300 timesteps, and plots the average virus population size as a
    function of time.

    numViruses: number of SimpleVirus to create for patient (an integer)
    maxPop: maximum virus population for patient (an integer)
    maxBirthProb: Maximum reproduction probability (a float between 0-1)        
    clearProb: Maximum clearance probability (a float between 0-1)
    numTrials: number of simulation runs to execute (an integer)
    """
    totalTime = 300
    noOfVirus = [0.0 for step in range(totalTime)]

    for trial in range(numTrials):
        viruses = [SimpleVirus(maxBirthProb, clearProb) for i in range(numViruses)]
        patient = Patient(viruses, maxPop)

        for step in range(totalTime):
            noOfVirus[step] += patient.update()

    for step in range(totalTime):
        noOfVirus[step] /= numTrials

    pylab.plot(range(totalTime), noOfVirus)
    pylab.title('Virus simulation without Drug')
    pylab.legend(['Virus without Drug'])
    pylab.xlabel('Time step')
    pylab.ylabel('Number of Viruses')
    pylab.show()

def plotFeaturePDF(ift, pft, outbase, fmin=0.0, fmax=1.0, fstep=0.01):
    """
    Plot a comparison between the input feature distribution and the 
    feature distribution of the predicted halos
    """
    plt.clf()
    nfbins = ( fmax - fmin ) / fstep
    fbins = np.logspace( fmin, fmax, nfbins )
    fcen = ( fbins[:-1] + fbins[1:] ) / 2

    plt.xscale( 'log', nonposx='clip' )
    plt.yscale( 'log', nonposy='clip' )
    
    ic, e, p = plt.hist( ift, fbins, label='Original Halos', alpha=0.5, normed=True )
    pc, e, p = plt.hist( pft, fbins, label='Added Halos', alpha=0.5, normed=True )

    plt.legend()
    plt.xlabel( r'$\delta$' )
    plt.savefig( outbase+'_fpdf.png' )

    fdtype = np.dtype( [ ('fcen', float), ('ifcounts', float), ('pfcounts', float) ] )
    fd = np.ndarray( len(fcen), dtype = fdtype )
    fd[ 'mcen' ] = fcen
    fd[ 'imcounts' ] = ic
    fd[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_fpdf.fit', fd )

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 plot_runtime_results(results):
    plt.rcParams["figure.figsize"] = 7,7
    plt.rcParams["font.size"] = 22
    matplotlib.rc("xtick", labelsize=24)
    matplotlib.rc("ytick", labelsize=24)

    params = {"text.fontsize" : 32,
              "font.size" : 32,
              "legend.fontsize" : 30,
              "axes.labelsize" : 32,
              "text.usetex" : False
              }
    plt.rcParams.update(params)
    
    #plt.semilogx(results[:,0], results[:,3], 'r-x', lw=3)
    #plt.semilogx(results[:,0], results[:,1], 'g-D', lw=3)
    #plt.semilogx(results[:,0], results[:,2], 'b-s', lw=3)

    plt.plot(results[:,0], results[:,3], 'r-x', lw=3, ms=10)
    plt.plot(results[:,0], results[:,1], 'g-D', lw=3, ms=10)
    plt.plot(results[:,0], results[:,2], 'b-s', lw=3, ms=10)

    plt.legend(["Chain", "Tree", "FFT Tree"], loc="upper left")
    plt.xticks([1e5, 2e5, 3e5])
    plt.yticks([0, 60, 120, 180])

    plt.xlabel("Problem Size")
    plt.ylabel("Runtime (sec)")
    return results

def fit_plot_unlabeled_data(unlabeled_data_x, labeled_data_x, labeled_data_y, fit_order, data_type, other_data_list, other_data_name):
    output = open('predictions.csv','wb')
    coeffs = np.polyfit(labeled_data_x, labeled_data_y, fit_order) #does poly git to nth deg on labeled data
    fit_eq = np.poly1d(coeffs) #Eqn from fit
    predicted_y = fit_eq(unlabeled_data_x)
    i = 0
    writer = csv.writer(output,delimiter=',')
    header = [str(data_type),str(other_data_name),'Predicted_Num_Inc']
    writer.writerow(header)
    while i < len(predicted_y):
        output_data = [unlabeled_data_x[i],other_data_list[i],predicted_y[i]]
        writer.writerow(output_data)
        print 'For '+str(data_type)+' of: '+str(unlabeled_data_x[i])+', Predicted Number of Incidents is: '+str(predicted_y[i])
        i = i + 1
    plt.scatter(unlabeled_data_x, predicted_y, color='blue', label='Predicted Number of Incidents')
    fit_line_x = np.arange(min(unlabeled_data_x), max(unlabeled_data_x), 1)
    plt.plot(fit_line_x, fit_eq(fit_line_x), color='red',linestyle='dashed',label=' Order '+str(fit_order)+' Polynomial Fit')
#____Use below line to plot actual data also!! 
    #plt.scatter(labeled_data_x, labeled_data_y, color='green', label='Actual Incident Report Data')
    plt.title('Predicted Number of 311 Incidents by '+str(data_type))
    plt.xlabel(str(data_type))
    plt.ylabel('Number of 311 Incidents')
    plt.grid()
    plt.xlim([min(unlabeled_data_x)-1500, max(unlabeled_data_x)+1500])
    plt.legend(loc='upper left')
    plt.show()

    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)