def plot_data(comp, c='b'):
    """utility function to make the Kantrowitz Limit Plot""" 

    MN = []
    W_tube = []
    W_kant = []

    for m in np.arange(.1,1.1,.1): 
        comp.Mach_pod = m
        comp.run()
        #print comp.radius_tube, comp.Mach_pod, comp.W_tube, comp.W_kant, comp.W_excess

        MN.append(m)
        W_kant.append(comp.W_kant)
        W_tube.append(comp.W_tube)

    fig = p.plot(MN,W_tube, '-', label="%3.1f Req."%(comp._tube_area/comp._inlet_area), lw=3, c=c)
    p.plot(MN,W_kant, '--', label="%3.1f Limit"%(comp._tube_area/comp._inlet_area),   lw=3, c=c)
    #p.legend(loc="best")
    p.tick_params(axis='both', which='major', labelsize=15)
    p.xlabel('Pod Mach Number', fontsize=18)
    p.ylabel('Flow Rate (kg/sec)', fontsize=18)
    p.title('Tube Flow Limits for Three Area Ratios', fontsize=20)

    return fig

    def plot(self, ylog10scale = False, timescale = "years", year = 25):
        """
        Generate figure and axis for the population structure
        timescale choose from "2N0", "4N0", "generation" or "years"
        """
        time = self.Time
        pop  = self.pop
        for i in range(1,len(self.pop)):
            if type(pop[i]) == type(""):
                # ignore migration commands, and replace by (unchanged) pop size
                pop[i] = pop[i-1]
        
        if time[0] != 0 :
            time.insert(0, float(0))
            pop.insert(0, float(1))
        
        if timescale == "years":
            time = [ti * 4 * self.scaling_N0 * year for ti in time ]
            pl.xlabel("Time (years, "+`year`+" years per generation)",  fontsize=20)    
            #pl.xlabel("Years")    
        elif timescale == "generation":
            time = [ti * 4 * self.scaling_N0 for ti in time ]
            pl.xlabel("Generations)")    
        elif timescale == "4N0":
            time = [ti*1 for ti in time ]
            pl.xlabel("Time (4N generations)")    
        elif timescale == "2N0":
            time = [ti*2 for ti in time ]
            pl.xlabel("Time (2N generations)")       
        else:
            print "timescale must be one of \"4N0\", \"generation\", or \"years\""
            return
        
        time[0] = time[1] / float(20)
        
        time.append(time[-1] * 2)
        yaxis_scaler = 10000
        
        pop = [popi * self.scaling_N0 / float(yaxis_scaler) for popi in pop ]
        pop.insert(0, pop[0])               
        pl.xscale ('log', basex = 10)        
        #pl.xlim(min(time), max(time))
        pl.xlim(1e3, 1e7)
        
        if ylog10scale:
            pl.ylim(0.06, 10000)
            pl.yscale ('log', basey = 10)            
        else:
            pl.ylim(0, max(pop)+2)
        
        pl.ylim(0,5)            
        pl.tick_params(labelsize=20)

        #pl.step(time, pop , color = "blue", linewidth=5.0)
        pl.step(time, pop , color = "red", linewidth=5.0)
        pl.grid()
        #pl.step(time, pop , color = "black", linewidth=5.0)
        #pl.title ( self.case + " population structure" )
        #pl.ylabel("Pop size ($*$ "+`yaxis_scaler` +")")
        pl.ylabel("Effective population size",fontsize=20 )

    def plot_currents(self, plot_which):

        global figures

        if plot_which == "all":
            plot_which = [i for i in range(1, self.instance["no_feeders"] + 1)]
        time_scale = [datetime.datetime(2000, 1, 1, 0) + datetime.timedelta(minutes=i) for i in range(1440)]
        plt.figure(figures)
        figures += 1
        lines_to_plot = len(plot_which)

        for i in range(lines_to_plot):

            CR = np.array(self.extract_from_csv(plot_which[i], self.instance["iteration"]))
            CR = np.reshape(CR, (-1, 3))

            plt.subplot(lines_to_plot, 1, 1 + i)
            if (i + 1) != lines_to_plot:
                plt.plot(time_scale, CR)
                plt.tick_params(axis='x', which='both', bottom='off', top='off', labelbottom='off')
            else:
                plt.plot(time_scale, CR)

        plt.subplot(lines_to_plot, 1, 1)
        plt.title("Current at the head of the Feeders")

def chart(SW, a, b, label, folder, FILE):
    pylab.ioff()
    fig_width_pt = 350 					     # Get this from LaTeX using \showthe\columnwidth
    inches_per_pt = 1.0/72.27                # Convert pt to inch
    golden_mean = ((5**0.5)-1.0)/2.0         # Aesthetic ratio
    fig_width = fig_width_pt*inches_per_pt   # width in inches   
    fig_height = fig_width*golden_mean       # height in inches
    fig_size =  [fig_width,fig_height]

    params = { 'backend': 'ps',
           'axes.labelsize': 10,
           'text.fontsize': 10,
           'legend.fontsize': 10,
           'xtick.labelsize': 8,
           'ytick.labelsize': 8,
           'text.usetex': True,
           'figure.figsize': fig_size }

    pylab.rcParams.update(params)

    home = '/home/nealbob'
    img_ext = '.pdf'

    pylab.figure()
    pylab.boxplot([SW['SWA'], SW['OA'], SW['NS']], whis=5)
    pylab.axhline(y=1.0, color='0.5', linewidth=0.5, alpha=0.75, linestyle=':')
    pylab.ylim(a, b)
    pylab.ylabel(label)
    pylab.tick_params(axis='x', which = 'both', labelbottom='off')
    pylab.figtext(0.225, 0.06, 'SWA', fontsize = 10)
    pylab.figtext(0.495, 0.06, 'OA', fontsize = 10)
    pylab.figtext(0.76, 0.06, 'NS', fontsize = 10)
    pylab.savefig(home + folder + FILE + img_ext)
    pylab.show()

def chart(idx, a, b, label, FILE):
    pylab.ioff()
    fig_width_pt = 350 					     # Get this from LaTeX using \showthe\columnwidth
    inches_per_pt = 1.0/72.27                # Convert pt to inch
    golden_mean = ((5**0.5)-1.0)/2.0         # Aesthetic ratio
    fig_width = fig_width_pt*inches_per_pt   # width in inches   
    fig_height = fig_width*golden_mean       # height in inches
    fig_size =  [fig_width*0.42,fig_height]

    params = { 'backend': 'ps',
           'axes.labelsize': 10,
           'text.fontsize': 10,
           'legend.fontsize': 10,
           'xtick.labelsize': 8,
           'ytick.labelsize': 8,
           'text.usetex': True,
           'figure.figsize': fig_size }

    pylab.rcParams.update(params)

    home = '/home/nealbob'
    folder = '/Dropbox/Thesis/IMG/chapter3/'
    img_ext = '.pdf'

    pylab.figure()
    pylab.boxplot(idx, whis=100)
    pylab.ylim(a, b)
    #pylab.ylabel(label)
    pylab.tick_params(axis='x', which = 'both', labelbottom='off')
    pylab.savefig(home + folder + FILE + img_ext)
    pylab.show()

def plots_of_lp_event_exchanges():
    pylab.title('Remote Events Sent Between LPs')
    data = np.loadtxt("analysisData/eventsExchanged-remote.csv", dtype=np.float_, delimiter = ",", skiprows=2, usecols=(2,3,4,5))
    outFile = outDir + 'countsOfLpToLpEventExchanges'
    pylab.plot(data[data[:,0].argsort()][:,0].astype(np.intc))
#    pylab.xlabel('Number of Events')
    pylab.tick_params(axis='x',labelbottom='off')
    pylab.ylabel('Number of Events Sent')
    display_graph(outFile)

    pylab.title('Timestamp Deltas of Remote Events')
    outFile = outDir + 'timeStampDeltasOfRemoteEvents'
    stride = max(int(max(len(data[:,1]),len(data[:,2]),len(data[:,3]))/20),1)
    pylab.plot(data[data[:,1].argsort()][:,1], color=colors[0], label="Minimum", marker='o', markevery=stride)
    pylab.plot(data[data[:,3].argsort()][:,3], color=colors[1], label="Average", marker='x', markevery=stride)
#    pylab.plot(data[data[:,2].argsort()][:,2], color=colors[2], label="Maximum", marker='*', markevery=stride)
    pylab.tick_params(axis='x',labelbottom='off')
    pylab.ylabel('Timestamp Delta (ReceiveTime - SendTime)')
    pylab.ylim([-.1,np.amax(data[:,3].astype(np.intc))+1])
#    pylab.yscale('log')
    pylab.legend(loc='best')
    display_graph(outFile)

    pylab.title('Histogram of Timestamp Deltas of Remote Events')
    outFile = outDir + 'timeStampDeltasOfRemoteEvents-hist'
    pylab.hist((data[:,1],data[:,3],data[:,2]), label=('Minimum', 'Average', 'Maximum'), color=(colors[0], colors[1], colors[2]), bins=10)
    pylab.xlabel('Timestamp Delta (ReceiveTime - SendTime)')    
    pylab.ylabel('Number of LPs')
    pylab.legend(loc='best')
    display_graph(outFile)

    return

def plotBar(data=None,color_id=None,figure_id=None,name=None,flag=False):
	ax = pl.subplot(figure_id)
	width = 0.8
	x=sp.arange(7)
	if not (name=="VaribenchSelected"):
		pl.bar(x-0.4,data,width=width,color=color_t[color_id],hatch="/o/o/")
	else:
		pl.bar(x-0.4,data,width=width,color=color_t[color_id],hatch="ooo")

	tmp = data.copy()
	tmp[1::] = 0
	pl.xticks(x,['All','Pure',']0.0,1.0[','[0.1,0.9]','[0.2,0.8]','[0.3,0.7]','[0.4,0.6]'],fontsize=font_size,rotation=90)
	ln = sp.log10(len(name))
	pl.text(3.5-ln,0.95,name)
	if flag:
		remove_border(left=False)
		pl.yticks([0.5,0.6,0.7,0.8,0.9,1.0])
		pl.grid(axis='y')
		pl.tick_params(axis='y',which="both",labelleft='off',left='off')
	else:
		pl.ylabel("AUC")
		remove_border()
		pl.yticks([0.5,0.6,0.7,0.8,0.9,1.0])
		pl.grid(axis='y')
	pl.ylim(0.5,1)
	pl.xlim(-0.5,7.5)
	return ax

def fixup_adjoint(current_data):
    import pylab
    size = 36
    addgauges(current_data)
    pylab.title('Adjoint Pressure', fontsize=size)
    pylab.xticks([-2, 0, 2, 4, 6], fontsize=size)
    pylab.tick_params(axis='y', labelleft='off')

def embedSystemtsne(embedding):
    pl.figure()
    pl.scatter(embedding[:,0],embedding[:,1],c=range(len(embedding[:,0])),linewidths=0)
    pl.title('t-distributed stochastic neighbor embedding of observed system states')
    pl.tick_params(labelleft='off', labelbottom='off')
    pl.colorbar().set_label('Time (ms)')
    pl.show()

def embedIndividualstsne(embedding):
    pl.figure()
    pl.scatter(embedding[:,0],embedding[:,1],c=d,linewidths=0)
    pl.title('t-distributed stochastic neighbor embedding of all individual neurons')
    pl.tick_params(labelleft='off', labelbottom='off')
    pl.colorbar().set_label('Parameter d (after-spike increment value of recovery variable u)')
    pl.show()

    def histplot(self, extradataA = [], extradataG = [], intensity = []):
        pylab.figure(figsize = (25,8))
        cat = ['NT, 500ng/mL DOX', 'DLG siRNA, 500ng/mL DOX', 'NuMA siRNA, 500ng/mL DOX', 'NT, 1ug/mL DOX']
        pops = []
        for i in xrange(3):
            pylab.subplot(1,3,i+1)
            pop = self.angles[(self.categories == i)]# &  (self.GFP > -np.log(12.5))]# & (intensity == 'r')]
            print "cat {0}, pop {1}, pop + GFP {2}".format(i, len(self.angles[self.categories == i]), len(pop))
            pops.append(pop)
            hist, binedges = np.histogram(pop, bins = 18)
            pylab.tick_params(axis='both', which='major', labelsize=25)
            pylab.plot(binedges[:-1], np.cumsum(hist)/1./len(pop), data.colors[i], label = data.cat[i], linewidth = 4)
            if len(extradataA) > i:
                print extradataA[i]
                h, bins = np.histogram(extradataA[i], bins= 18)
                hbis = h/1./len(extradataA[i])
                x, y = [], []
                for index in xrange(len(hbis)):
                    x.extend([bins[index], bins[index+1]])
                    y.extend([hbis[index], hbis[index]])
                print x, y, len(x)
                pylab.tick_params(axis='both', which='major', labelsize=25)
                pylab.plot(bins[:-1], np.cumsum(h)/1./len(extradataA[i]), 'k', linewidth = 4)

            pylab.xlabel("Angle (degre)", fontsize = 25)
            #pylab.title(cat[i])
            pylab.ylim([0., 1.2])
            pylab.legend(loc = 2, prop = {'size' : 20})
        for ip, p in enumerate(pops):
            for ip2, p2 in enumerate(pops):
                ksstat, kspval = scipy.stats.ks_2samp(p2, p)
                print "#### cat{0} & cat{3} : ks Stat {1}, pvalue {2}".format(ip, ksstat, kspval, ip2)
        pylab.show()

  def make_plotII(self):
    # retrieve data
    D=self.D

    kmap={}
    kmap['Q2 = 2']   = {'c':'r','ls':'-'}
    kmap['Q2 = 5']   = {'c':'g','ls':'--'}
    kmap['Q2 = 10']  = {'c':'b','ls':'-.'}
    kmap['Q2 = 100'] = {'c':'k','ls':':'}


    ax=py.subplot(111)
    DF=D['AV18']
    for Q2 in [2,5,10,100]:
      k='Q2 = %d'%Q2
      Q2=float(k.split('=')[1])
      DF=D['AV18'][D['AV18'].Q2==Q2]
      cls=kmap[k]['c']+kmap[k]['ls']
      ax.plot(DF.X,DF.THEORY,cls,lw=2.0,label=r'$Q^2=%0.0f~{\rm GeV}^2$'%Q2)

    ax.set_xlabel('$x$',size=25)
    ax.set_ylabel(r'$F_2^d\, /\, F_2^N$',size=25)
    ax.set_ylim(0.97,1.08)
    ax.axhline(1,color='k',ls='-',alpha=0.2)

    ax.legend(frameon=0,loc=2,fontsize=22)
    py.tick_params(axis='both',labelsize=22)
    py.tight_layout()
    py.savefig('gallery/F2d_F2_II.pdf')

  def make_plotI(self):
    # retrieve data
    D=self.D

    kmap={}
    kmap['AV18']   = {'c':'r','ls':'-'}
    kmap['CDBONN'] = {'c':'g','ls':'--'}
    kmap['WJC1']   = {'c':'k','ls':'-.'}
    kmap['WJC2']   = {'c':'b','ls':':'}

    ax=py.subplot(111)
    for k in ['AV18','CDBONN','WJC1','WJC2']:
      DF=D[k]
      DF=DF[DF.Q2==10]
      if k=='CDBONN':
        label='CDBonn'
      else:
        label=k
      cls=kmap[k]['c']+kmap[k]['ls']
      ax.plot(DF.X,DF.THEORY,cls,lw=2.0,label=tex(label))

    ax.set_xlabel('$x$',size=25)
    ax.set_ylabel(r'$F_2^d\, /\, F_2^N$',size=25)
    ax.set_ylim(0.97,1.08)
    ax.axhline(1,color='k',ls='-',alpha=0.2)

    ax.legend(frameon=0,loc=2,fontsize=22)
    py.tick_params(axis='both',labelsize=22)
    py.tight_layout()
    py.savefig('gallery/F2d_F2_I.pdf')
    py.close()

    def createCoefGraph(data, nFig, lim, ymin):
        plt.figure(nFig)
        plt.suptitle('Coef')
        nBase = data.shape[0]
        nSubCols = nBase / 10
        if nSubCols > 0:
            nSubRows = nBase / nSubCols
        else:
            nSubRows = nBase 
            nSubCols = 	1
        # print data.shape

        # サンプリング周波数とシフト幅によって式を変える必要あり
        timeLine = [i * 1024 / 8000.0 for i in range(data.shape[1])]
        # print len(timeLine)
        for i in range(nBase):
            plt.subplot(nSubRows, nSubCols, i + 1)
            plt.tick_params(labelleft='off', labelbottom='off')
            # FIXME: Arguments of X
            # plt.plot(timeLine, data[i,:])
            if lim:
                plt.ylim(ymin=ymin)
            plt.plot(timeLine, data[i,:])
        # Beacuse I want to add lable in bottom, xlabel is declaration after loop.
        plt.tick_params(labelleft='off', labelbottom="on")
        plt.xlabel('time [ms]')

def profile_of_local_events_exec_by_lp():
    pylab.title('Locally Generated Events')
    outFile = outDir + 'percentOfExecutedEventsThatAreLocal'
    data = np.loadtxt("analysisData/eventsExecutedByLP.csv", dtype=np.intc, delimiter = ",", skiprows=2, usecols=(1,2,3))
    x_index = np.arange(len(data))
    pylab.plot(x_index, sorted(percent_of_LP_events_that_are_local(data)))
    pylab.xlabel('LPs (sorted by percent local)')
    pylab.tick_params(axis='x',labelbottom='off')
    pylab.ylabel('Percent of Total Executed (Ave=%.2f%%)' % np.mean(percent_of_LP_events_that_are_local(data)))
    pylab.ylim((0,100))
    # fill the area below the line
    ax = pylab.gca()
#    ax.fill_between(x_index, sorted(percent_of_LP_events_that_are_local(data)), 0, facecolor=colors[0])
    ax.get_yaxis().set_major_formatter(mpl.ticker.FormatStrFormatter('%.1f%%'))
    display_graph(outFile)

    pylab.title('Locally Generated Events Executed')
    outFile = outDir + 'percentOfExecutedEventsThatAreLocal-histogram'
    pylab.hist(sorted(percent_of_LP_events_that_are_local(data)))
    ax = pylab.gca()
    ax.get_xaxis().set_major_formatter(mpl.ticker.FormatStrFormatter('%.1f%%'))
    pylab.xlabel('Percent of Local Events Executed')
    pylab.ylabel('Number of LPs (Total=%s)' % "{:,}".format(total_lps))
    display_graph(outFile)
    return

def plot_data(p, c='b'):
    '''utility function to make the Kantrowitz Limit Plot'''
    Machs = []
    W_tube = []
    W_kant = []
    for Mach in np.arange(.2, 1.1, .1):
        p['comp.Mach'] = Mach
        p.run()
        Machs.append(Mach)
        W_kant.append(p['comp.W_kant'])
        W_tube.append(p['comp.W_tube'])
    print('Area in:', p['comp.inlet.area_out'])
    fig = pylab.plot(Machs,
                     W_tube,
                     '-',
                     label="%3.1f Req." %
                     (p['comp.tube_area'] / p['comp.inlet.area_out']),
                     lw=3,
                     c=c)
    pylab.plot(Machs,
               W_kant,
               '--',
               label="%3.1f Limit" %
               (p['comp.tube_area'] / p['comp.inlet.area_out']),
               lw=3,
               c=c)
    pylab.tick_params(axis='both', which='major', labelsize=15)
    pylab.xlabel('Pod Mach Number', fontsize=18)
    pylab.ylabel('Flow Rate (kg/sec)', fontsize=18)
    pylab.title('Tube Flow Limits for Three Area Ratios', fontsize=20)
    return fig

    def createBasisGraph(self, data):
        plt.figure(self.nFig)
        plt.suptitle('Basis')
        nBase = data.shape[1]
        # The cols number of subplot is 
        nSubCols = nBase / 10
        if nSubCols > 0:
            if nBase % 2 == 0:
                nSubRows = nBase / nSubCols
            else:
                nSubRows = nBase / nSubCols + 1
        else:
            nSubRows = nBase 
            nSubCols = 1
        # freqList = np.fft.fftfreq(513, d = 1.0 / 44100)

        for i in range(nBase):
            nowFig = self.nFig + (i / nSubRows) + 1
            # Because Index of graph is start by 1, The Graph index start from i + 1.
            plt.subplot(nSubRows, nSubCols, i + 1)
            plt.tick_params(labelleft='off', labelbottom='off')

            # FIXME
            #plt.ylabel(self.st5[i%12] + str(i/12  + 1))
            plt.ylabel(str(i))
            plt.plot(data[:,i])
        # Beacuse I want to add lable in bottom, xlabel is declaration after loop.
        plt.tick_params(labelleft='off', labelbottom='on')
        plt.xlabel('frequency [Hz]')

        #self.nFig += nowFig
        self.nFig += 1 

    def plotAgainstGFP(self, extradataA = [], extradataG = [], intensity = [], seq = []):
        fig1 = pylab.figure(figsize = (25, 10))
        print len(self.GFP)
        for i in xrange(min(len(data.cat), 3)):
            print len(self.GFP[self.categories == i])
            vect = []
            pylab.subplot(1,3,i+1)
            #pylab.hist(self.GFP[self.categories == i], bins = 20, color = data.colors[i])
            pop = self.GFP[self.categories == i]
            pylab.plot(self.GFP[self.categories == i], self.angles[self.categories == i], data.colors[i]+'o', markersize = 8)#, label = data.cat[i])
            print "cat", i, "n pop", len(self.GFP[(self.categories == i) & (self.GFP > -np.log(12.5))])
            x = np.linspace(np.min(self.GFP[self.categories == i]), np.percentile(self.GFP[self.categories == i], 80),40)
            #fig1.canvas.mpl_connect('pick_event', onpick)
            for j in x:
                vect.append(np.median(self.angles[(self.GFP > j) & (self.categories == i)]))

            pylab.plot([-4.5, -0.5], [vect[0], vect[0]], data.colors[i], label = "mediane de la population entiere", linewidth = 5)
            print vect[0], vect[np.argmax(x > -np.log(12.5))]
            pylab.plot([-np.log(12.5), -0.5], [vect[np.argmax(x > -np.log(12.5))] for k in  [0,1]], data.colors[i], label = "mediane de la population de droite", linewidth = 5, ls = '--')
            pylab.axvline(x = -np.log(12.5), color = 'm', ls = '--', linewidth = 3)
            pylab.xlim([-4.5, -0.5])
            pylab.legend(loc = 2, prop = {'size':17})

            pylab.title(data.cat[i].split(',')[0], fontsize = 24)
            pylab.xlabel('score GFP', fontsize = 20)
            pylab.ylabel('Angle (degre)', fontsize = 20)
            pylab.tick_params(axis='both', which='major', labelsize=20)
            pylab.ylim([-5, 105])
            ##pylab.xscale('log')
        pylab.show()

 def GraphMaker(self,name,c,d,e):
     """Produces graphs based on the user entered equations in the Result Window"""
     
     model = name.rstrip("']")
     model = model.lstrip("'[")
     model = str(model)
     name = ''
     for i in xrange(len(model)):
         if model[i] == '_':
             name += ' '
         else:
             name += model[i]
             
     __location__ = os.path.dirname(sys.argv[0]) 
     loc = os.path.join(__location__, 'Data/')
     
     Params = ["w_0", "w_a", "w_p", "w_DE", "Omega_M", "Omega_DE", "x1", "x2", "y1", "y2", "c1", "c2", "u", "Lambda_1", "Lambda_2", "n"]
     for entry in xrange(len(Params)):
         exec( Params[entry] + " = np.genfromtxt(loc+str(model)+'.dat', usecols = "+str(entry)+", skip_header = 1)" )
             
     A = eval(c)
     B = eval(d)
         
     py.figure(int(e))
     py.plot(A, B, 'b.')
     py.xlabel(c, fontsize = 55)
     py.ylabel(d, fontsize = 55)
     py.title('Results from '+str(name), fontsize = 55)
     py.tick_params(labelsize = 35, size = 15, width = 5, top = 0, right = 0)
         
     py.show()

def fixup_innerprod(current_data):
    import pylab
    size = 28
    addgauges(current_data)
    pylab.title('Inner Product', fontsize=size)
    pylab.xticks(fontsize=size)
    pylab.tick_params(axis='y', labelleft='off')