def plot_single_output_multiple_days(self, single_output, days_list=SOLSTICES, title_str=''):
        # Output: graph with subplots for each day
        # Inputs: df_dict: dict of dataframes of PVsyst result
        #        single_output: keyword name of data to plot
        #        days_list: list of datetimes as strings in format %Y-%m-%d; used for subplots' titles
        single_output_df = self.get_single_output_df(single_output)
        legend_list = list(single_output_df.columns)
        y_info = get_y_axis_info(single_output)
        plt.figure(figsize=(8, 15))
        plt.suptitle(self.location + ': ' + single_output + '\n' + title_str, fontsize=16)

        ax = []
        for d in range(1, len(days_list) + 1):
            if d == 1:
                ax.append(plt.subplot(len(days_list), 1, d))
                plt.ylabel(y_info[UNITS], fontsize=14)
            else:
                ax.append(plt.subplot(len(days_list), 1, d, sharex=ax[0], sharey=ax[0]))
            plt.grid(True, which='major', axis='x')
            ax[d - 1].plot(single_output_df[days_list[d - 1]].index.hour, single_output_df[days_list[d - 1]]/y_info[SCALE],
                           marker='.')
            # if show_flat == True:
            #     ax[d - 1].plot(flat_df[days_list[d - 1]].index.hour, flat_df[days_list[d - 1]], 'k--')
            ax[d - 1].legend(legend_list, loc='lower right', fontsize=9)
            ax[d - 1].set_title(days_list[d - 1], fontsize=12)

        plt.xlabel('Hour of day', fontsize=14)
        ax[0].set_xlim([0, 24])
        ax[0].set_xticks([0, 6, 12, 18, 24])
        ax[0].set_ylim(y_info[LIM])

	def plotPSD(self, chan, time_interval):
		Npackets = np.int(time_interval * self.accum_freq)
		plot_range = (Npackets / 2) + 1
		figure = plt.figure(num= None, figsize=(12,12), dpi=80, facecolor='w', edgecolor='w')
		# I 
		plt.suptitle('Channel ' + str(chan) + ' , Freq = ' + str((self.freqs[chan] + self.LO_freq)/1.0e6) + ' MHz') 
		plot1 = figure.add_subplot(311)
		plot1.set_xscale('log')
		plot1.set_autoscale_on(True)
		plt.ylim((-160,-80))
		plt.title('I')
		line1, = plot1.plot(np.linspace(0, self.accum_freq/2., (Npackets/2) + 1), np.zeros(plot_range), label = 'I', color = 'green', linewidth = 1)
		plt.grid()
		# Q
		plot2 = figure.add_subplot(312)
		plot2.set_xscale('log')
		plot2.set_autoscale_on(True)
		plt.ylim((-160,-80))
		plt.title('Q')
		line2, = plot2.plot(np.linspace(0, self.accum_freq/2., (Npackets/2) + 1), np.zeros(plot_range), label = 'Q', color = 'red', linewidth = 1)
		plt.grid()
		# Phase
		plot3 = figure.add_subplot(313)
		plot3.set_xscale('log')
		plot3.set_autoscale_on(True)
		plt.ylim((-120,-70))
		#plt.xlim((0.0001, self.accum_freq/2.))
		plt.title('Phase')
		plt.ylabel('dBc rad^2/Hz')
		plt.xlabel('log Hz')
		line3, = plot3.plot(np.linspace(0, self.accum_freq/2., (Npackets/2) + 1), np.zeros(plot_range), label = 'Phase', color = 'black', linewidth = 1)
		plt.grid()
		plt.show(block = False)
		count = 0
		stop = 1.0e10
		while count < stop:
			Is, Qs, phases = self.get_stream(chan, time_interval)
			I_mags = np.fft.rfft(Is, Npackets)
			Q_mags = np.fft.rfft(Is, Npackets)
			phase_mags = np.fft.rfft(phases, Npackets)
			I_vals = (np.abs(I_mags)**2 * ((1./self.accum_freq)**2 / (1.0*time_interval)))
			Q_vals = (np.abs(Q_mags)**2 * ((1./self.accum_freq)**2 / (1.0*time_interval)))
			phase_vals = (np.abs(phase_mags)**2 * ((1./self.accum_freq)**2 / (1.0*time_interval)))
			phase_vals = 10*np.log10(phase_vals)
			phase_vals -= phase_vals[0]
			#line1.set_ydata(Is)
			#line2.set_ydata(Qs)
			#line3.set_ydata(phases)
			line1.set_ydata(10*np.log10(I_vals))
			line2.set_ydata(10*np.log10(Q_vals))
			line3.set_ydata(phase_vals)
			plot1.relim()
			plot1.autoscale_view(True,True,False)
			plot2.relim()
			plot2.autoscale_view(True,True,False)
			#plot3.relim()
			plot3.autoscale_view(True,True,False)
			plt.draw()
			count +=1
		return

def scatter_time_vs_s(time, norm, point_labels, title):
    plt.figure()
    size = 100
    for i, l in enumerate(sorted(norm.keys())):
        if l is not "fbpca":
            plt.scatter(time[l], norm[l], label=l, marker='o', c='b', s=size)
            for label, x, y in zip(point_labels, list(time[l]), list(norm[l])):
                plt.annotate(label, xy=(x, y), xytext=(0, -80),
                             textcoords='offset points', ha='right',
                             arrowprops=dict(arrowstyle="->",
                                             connectionstyle="arc3"),
                             va='bottom', size=11, rotation=90)
        else:
            plt.scatter(time[l], norm[l], label=l, marker='^', c='red', s=size)
            for label, x, y in zip(point_labels, list(time[l]), list(norm[l])):
                plt.annotate(label, xy=(x, y), xytext=(0, 30),
                             textcoords='offset points', ha='right',
                             arrowprops=dict(arrowstyle="->",
                                             connectionstyle="arc3"),
                             va='bottom', size=11, rotation=90)

    plt.legend(loc="best")
    plt.suptitle(title)
    plt.ylabel("norm discrepancy")
    plt.xlabel("running time [s]")

    def display_channel_efficiency(self):

        size = 0

        start_time = self.pcap_file[0].time
        end_time = self.pcap_file[len(self.pcap_file) - 1].time

        duration = (end_time - start_time)/1000

        for i in range(len(self.pcap_file) - 1):
            size += len(self.pcap_file[i])
        ans = (((size * 8) / duration) / BW_STANDARD_WIFI) * 100
        ans = float("%.2f" % ans)
        labels = ['utilized', 'unutilized']
        sizes = [ans, 100.0 - ans]
        colors = ['g', 'r']

        # Make a pie graph
        plt.clf()
        plt.figure(num=1, figsize=(8, 6))
        plt.axes(aspect=1)
        plt.suptitle('Channel efficiency', fontsize=14, fontweight='bold')
        plt.title("Bits/s: " + str(float("%.2f" % ((size*8)/duration))),fontsize = 12)
        plt.rcParams.update({'font.size': 17})
        plt.pie(sizes, labels=labels, autopct='%.2f%%', startangle=60, colors=colors, pctdistance=0.7, labeldistance=1.2)

        plt.show()

    def display_graph_by_specific_mac(self, mac_address):

        G = nx.Graph()

        count = 0
        edges = set()
        edges_list = []


        for pkt in self.pcap_file:

            src = pkt[Dot11].addr1
            dst = pkt[Dot11].addr2

            if mac_address in [src, dst]:
                edges_list.append((src, dst))
                edges.add(src)
                edges.add(dst)

        plt.clf()
        plt.suptitle('Communicating with ' + str(mac_address), fontsize=14, fontweight='bold')
        plt.title("\n Number of Communicating Users: " + str(int(len(edges))))
        plt.rcParams.update({'font.size': 10})
        G.add_edges_from(edges_list)
        nx.draw(G, with_labels=True, node_color=MY_COLORS)
        plt.show()

def run_div_test(fld, exact, title='', show=False, ignore_inexact=False):
    t0 = time()
    result_numexpr = viscid.div(fld, preferred="numexpr", only=False)
    t1 = time()
    logger.info("numexpr magnitude runtime: %g", t1 - t0)

    result_diff = viscid.diff(result_numexpr, exact)['x=1:-1, y=1:-1, z=1:-1']
    if not ignore_inexact and not (result_diff.data < 5e-5).all():
        logger.warning("numexpr result is far from the exact result")
    logger.info("min/max(abs(numexpr - exact)): %g / %g",
                np.min(result_diff.data), np.max(result_diff.data))

    planes = ["y=0j", "z=0j"]
    nrows = 2
    ncols = len(planes)
    _, axes = plt.subplots(nrows, ncols, squeeze=False)

    for i, p in enumerate(planes):
        vlt.plot(result_numexpr, p, ax=axes[0, i], show=False)
        vlt.plot(result_diff, p, ax=axes[1, i], show=False)

    plt.suptitle(title)
    vlt.auto_adjust_subplots(subplot_params=dict(top=0.9))

    plt.savefig(next_plot_fname(__file__))
    if show:
        vlt.mplshow()

def template_matching():
    img = cv2.imread('messi.jpg',0)
    img2 = img.copy()
    template = cv2.imread('face.png',0)
    w, h = template.shape[::-1]

    # All the 6 methods for comparison in a list
    methods = ['cv2.TM_CCOEFF', 'cv2.TM_CCOEFF_NORMED', 'cv2.TM_CCORR',
            'cv2.TM_CCORR_NORMED', 'cv2.TM_SQDIFF', 'cv2.TM_SQDIFF_NORMED']

    for meth in methods:
        img = img2.copy()
        method = eval(meth)

        # Apply template Matching
        res = cv2.matchTemplate(img,template,method)
        min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res)

        # If the method is TM_SQDIFF or TM_SQDIFF_NORMED, take minimum
        if method in [cv2.TM_SQDIFF, cv2.TM_SQDIFF_NORMED]:
            top_left = min_loc
        else:
            top_left = max_loc
        bottom_right = (top_left[0] + w, top_left[1] + h)

        cv2.rectangle(img,top_left, bottom_right, 255, 2)

        plt.subplot(121),plt.imshow(res,cmap = 'gray')
        plt.title('Matching Result'), plt.xticks([]), plt.yticks([])
        plt.subplot(122),plt.imshow(img,cmap = 'gray')
        plt.title('Detected Point'), plt.xticks([]), plt.yticks([])
        plt.suptitle(meth)

        plt.show()

    def generate_line_chart(self, data, node_type, field, append_table=True):
        output = []
        common.bash("mkdir -p ../visualizer/include/pic")
        common.bash("mkdir -p ../visualizer/include/csv")
        common.printout("LOG","generate %s line chart" % node_type)
        for field_column, field_data in data.items():
            pyplot.figure(figsize=(9, 4))
            for node, node_data in field_data.items():
                pyplot.plot(node_data, label=node)
            pyplot.xlabel("time(sec)")
            pyplot.ylabel("%s" % field_column)
            # Shrink current axis's height by 10% on the bottom
            pyplot.legend(loc = 'center left', bbox_to_anchor = (1, 0.5), prop={'size':6})
            pyplot.grid(True)
            pyplot.suptitle("%s" % field_column)
            pic_name = '%s_%s_%s.png' % (node_type, field, re.sub('[/%]','',field_column))
            pyplot.savefig('../visualizer/include/pic/%s' % pic_name) 
            pyplot.close()
            line_table = []
            csv = self.generate_csv_from_json(field_data,'line_table',field_column)
            csv_name = '%s_%s_%s.csv' % (node_type, field, re.sub('[/%]','',field_column))
            with open( '../visualizer/include/csv/%s' % csv_name, 'w' ) as f:
                f.write( csv )
#output.append("<div class='cetune_pic' id='%s_%s_pic'><button><a href='./include/csv/%s'>Download detail csv table</a></button><img src='./include/pic/%s' alt='%s' style='height:400px; width:1000px'></div>" % (field, re.sub('[/%]','',field_column), csv_name, pic_name, field_column))
            output.append("<div class='cetune_pic' id='%s_%s_pic'><img src='./include/pic/%s' alt='%s' style='height:400px; width:1000px'><button><a href='./include/csv/%s'>Download detail csv table</a></button></div>" % (field, re.sub('[/%]','',field_column), pic_name, field_column, csv_name))
        return output

def run_mag_test(fld, title="", show=False):
    vx, vy, vz = fld.component_views()  # pylint: disable=W0612
    vx, vy, vz = fld.component_fields()

    try:
        t0 = time()
        mag_ne = viscid.magnitude(fld, preferred="numexpr", only=False)
        t1 = time()
        logger.info("numexpr mag runtime: %g", t1 - t0)
    except viscid.verror.BackendNotFound:
        xfail("Numexpr is not installed")

    planes = ["z=0", "y=0"]
    nrows = 4
    ncols = len(planes)

    _, axes = plt.subplots(nrows, ncols, sharex=True, sharey=True, squeeze=False)

    for ind, p in enumerate(planes):
        vlt.plot(vx, p, ax=axes[0, ind], show=False)
        vlt.plot(vy, p, ax=axes[1, ind], show=False)
        vlt.plot(vz, p, ax=axes[2, ind], show=False)
        vlt.plot(mag_ne, p, ax=axes[3, ind], show=False)

    plt.suptitle(title)
    vlt.auto_adjust_subplots(subplot_params=dict(top=0.9, right=0.9))
    plt.gcf().set_size_inches(6, 7)

    plt.savefig(next_plot_fname(__file__))
    if show:
        vlt.mplshow()

def localisationFault():
    # read file
    faults = []
    for line in f1:
        line = line.split('\t')
        lst = [float(x) for x in line]
        faults.append(lst)
    
    # extract
    xFault = list(map(lambda l : l[0], faults))
    yFault = list(map(lambda l : l[1], faults))
    thetaFault = list(map(lambda l : l[2], faults))
    
    # number of items
    num = len(xFault) + 1
    time = list(range(1, num))
    
    plt.suptitle('Robot localisation fault', fontsize=20)
    # plot
    plt.plot(time, xFault, 'r', label='x fault')
    plt.plot(time, yFault, 'b', label='y fault')
    plt.plot(time, thetaFault, 'y', label='theta fault')
    plt.legend(numpoints=3)
    # xMin - xMax and yMin - yMax
    plt.axis([0, num, 0, 1.2])
    plt.show()

def boxPositions():
    positons = []
    f2 = open(Stats.BOX_POSITIONS_FILE)
    for line in f2:
        line = line.split('\t')
        lst = [float(x) for x in line]
        positons.append(lst)
        
    approX = list(map(lambda l : l[0], positons))
    approY = list(map(lambda l : l[1], positons))
    trueX = list(map(lambda l : l[2], positons))
    trueY = list(map(lambda l : l[3], positons))
    
    # number of items
    num = len(approX) + 1
    time = list(range(1, num))
    
    plt.suptitle('Box positions', fontsize=20)
    # plot
    plt.subplot(211)
    plt.ylabel('x', fontsize=18)
    plt.plot(time, approX, 'r', label='Approximate Position')
    plt.plot(time, trueX, 'g', label='True Position')
    plt.legend(numpoints=1)
    plt.subplot(212)
    plt.ylabel('y', fontsize=18)
    plt.plot(time, approY, 'r')
    plt.plot(time, trueY, 'g')

    plt.show()

 def _drawTitle(self,isContour=False):
     #Add a title
     hlon = self.shakemap.getEventDict()['lon']
     hlat = self.shakemap.getEventDict()['lat']
     edict = self.shakemap.getEventDict()
     eloc = edict['event_description']
     timestr = edict['event_timestamp'].strftime('%b %d, %Y %H:%M:%S')
     mag = edict['magnitude']
     if hlon < 0:
         lonstr = 'W%.2f' % np.abs(hlon)
     else:
         lonstr = 'E%.2f' % hlon
     if hlat < 0:
         latstr = 'S%.2f' % np.abs(hlat)
     else:
         latstr = 'N%.2f' % hlat
     dep = edict['depth']
     eid = edict['event_id']
     net = edict['event_network']
     if not eid.startswith(net):
         eid = net + eid
     tpl = (timestr,mag,latstr,lonstr,dep,eid)
     layername = 'MMI'
     if isContour:
         layername = self.contour_layer.upper()
     plt.suptitle('USGS ShakeMap (%s): %s' % (layername,eloc),fontsize=14,verticalalignment='top',y=0.95)
     plt.title('%s UTC M%.1f %s %s Depth: %.1fkm ID:%s' % tpl,fontsize=10,verticalalignment='bottom')
     return eid

def plot_marginals(state_space,p,D,name,rank,t,to_file = False):
	import matplotlib
	#matplotlib.use("PDF")
	#matplotlib.rcParams['figure.figsize'] = 5,10
	import matplotlib.pyplot as pl
	pl.suptitle("time: "+ str(t)+" units")
	print("time : "+ str(t))
	for i in range(D):
		marg_X = np.unique(state_space[:,i])
		A = np.where(marg_X[:,np.newaxis] == state_space[:,i].T[np.newaxis,:],1,0)
		marg_p = np.dot(A,p)
		pl.subplot(int(D/2)+1,2,i+1)
		pl.plot(marg_X,marg_p)
		if to_file == True:
			string_state = str(t) 
			string_prob = str(t) 
			f = open("Visuals/cummulative_state_"+ str(i)+".txt",'a')
			g = open("Visuals/cummulative_prob_"+ str(i)+".txt",'a')
			for j in range(len(marg_X)):
				string_state +=','+ str(marg_X[j])
				string_prob += ','+ str(marg_p[j])
			string_state += '\n'
			string_prob += '\n'
			f.write(string_state)
			g.write(string_prob)
			f.close()
			g.close()	
	 

	pl.savefig("Visuals/marginal_"+name+".pdf",format='pdf')
	pl.show()
	pl.clf()

def compare_single_output_across_batches(batches_list, variants_list, params_list, single_output, days_list=SOLSTICES):
    # ORDER MATTERS for batches_list, variants_list, and params_list!!
    # e.g. this method will compare the param[i] of variant[i] of the batch[i] with
    # the param[i+1] of variant[i+1] of batch[i+1], etc.
    results_df_list = []
    legend_list = []
    for b in range(len(batches_list)):
        results_df_list.append(batches_list[b].results_dict.get(variants_list[b]))
        # legend_list.append(batches_list[b].location[:11] +
        #                    str(batches_list[b].variant_output_map[variants_list[b]]) +
        #                    batches_list[b].output_name)

    plt.figure(figsize=(8, 15))
    plt.suptitle(single_output, fontsize=16)
    y_info = get_y_axis_info(single_output)

    ax = []
    for d in range(1, len(days_list) + 1):
        if d == 1:
            ax.append(plt.subplot(len(days_list), 1, d))
            plt.ylabel(y_info[UNITS], fontsize=14)
        else:
            ax.append(plt.subplot(len(days_list), 1, d, sharex=ax[0], sharey=ax[0]))
        plt.grid(True, which='major', axis='x')
        [ax[d - 1].plot(r[days_list[d - 1]].index.hour, r.loc[days_list[d - 1], single_output]/y_info[SCALE])
         for r in results_df_list]
        ax[d - 1].set_title(days_list[d - 1], fontsize=12)

    plt.xlabel('Hour of day', fontsize=14)
    ax[0].set_xlim([0, 24])
    ax[0].set_xticks([0, 6, 12, 18, 24])
    ax[0].set_ylim(y_info[LIM])
    ax[0].legend(params_list, loc='lower right', fontsize=9)

def boxplotThem(dataToPlot,title):
    means = [np.mean(item) for item in dataToPlot]
    fig = plt.figure(1, figsize=(9,6))
    ax = fig.add_subplot(111)

    plt.rc('text', usetex=True)
    plt.rc('font', family='serif')
    for i in range(25):
        print('length ',str(i+1))
        for item in dataToPlot[i]:
            ax.scatter(i+1,item)

    ax.plot(list(range(1,26)),means,linewidth=4,color='r')

    ax.set_xticks([1,5,10,15,20,25])
    ax.set_xticklabels([r"$1",r"$5$",r"$10$",r"$15$",r"$20$",r"$25$"],fontsize = 25)

    ax.set_yscale('log')
    ax.set_yticks([0.00000001,0.000001,0.0001,0.01,1])
    ax.xaxis.set_tick_params(width=1.5)
    ax.yaxis.set_tick_params(width=1.5)
    ax.set_yticklabels([r"$10^{-8}$",r"$10^{-6}$",r"$10^{-4}$",r"$10^{-2}$",r"$1$"],fontsize = 25)
    ax.get_yaxis().get_major_formatter().labelOnlyBase = False
    ax.set_ylabel('relative population',fontsize = 30)
    ax.set_xlabel(r'length',fontsize = 30)
    ax.set_xlim(0,26)
    ax.set_ylim(0.00000005)
    plt.suptitle(title,fontsize=25)
    plt.savefig(r.outputDir+'distr.png')

def plot_TP():
    with open ('p_files/ROC_table_' + str(seg) + '_pc' + str(p) + '_k' + str(k) + '.p', 'rb') as f:
        ROC_table = pickle.load(f)
    with open ('p_files/species_stats.p', 'rb') as f:
        species_table = pickle.load(f)

    with open ('p_files/species.p', 'rb') as f:
        species_name = pickle.load(f)

    xes = []#[item['number'] for item in ROC_table.values()]
    yes = []#[item['fp'] for item in ROC_table.values()]
    label = []
    low_TP = []
    for specie in ROC_table:
        xes.append(species_table[specie])
        yes.append(ROC_table[specie]['tp_rate'])
        label.append(specie)
        if float(ROC_table[specie]['tp_rate']) < 0.3 and species_table[specie] > 100:
            #print(ROC_table[specie]['tp_rate'])
            low_TP.append((specie, ROC_table[specie]['tp']))

    fig, ax = plt.subplots()
    plt.subplots_adjust(bottom=0.1)
    ax.plot([0,max(xes)],[0.3,0.3], ls="--")
    ax.scatter(xes, yes, marker = '.')
    for i, txt in enumerate(label):
        if txt in interesting:
        #if float(yes[i]) < 0.3 and xes[i] > 100 :
            #print(txt)
            #ax.annotate(parser.get_specie_name('../data/train/', str(species_name[txt][0]) + '.xml'), (xes[i],yes[i]))
            ax.annotate(txt, (xes[i],yes[i]))
    plt.suptitle('False Positive', fontsize = 14)
#    ax.set_xlabel('Principal Components')
#    ax.set_ylabel('Percentage of Variance')
    plt.show()

def plotInLen(r,freqs,fits,title):
    nc = 3
    minLength=8
    maxLength=25
    fig, axes = plt.subplots(
        nrows=int((maxLength-minLength+1)/2/nc), ncols=nc, figsize=(12, 12)
        )
    index = 0
    for (length,freqs) in freqs.items():
        if length>=minLength and length<=maxLength and length%2==1:
            sortedFreqs = OrderedDict(sorted(freqs.items(),key=lambda t: t[0], reverse=False))
            axes[int((index)/nc),(index)%(nc)].plot(
                list(sortedFreqs.keys()),list(sortedFreqs.values()),'o'
                )
            m, c = fits[length-1]
            axes[int((index)/nc),(index)%(nc)].plot(
                list(sortedFreqs.keys()),
                [(10**(c)*xi**(m)) for xi in sortedFreqs.keys()],
                linewidth = 3,
                label = 'y = '+str("%.2f" %(10**c))+'x^'+str("%.2f" %m)
                )
            axes[int((index)/nc),(index)%(nc)].set_yscale('log')
            axes[int((index)/nc),(index)%(nc)].set_xscale('log')
            #axes[int((index)/nc),(index)%(nc)].legend()
            axes[int((index)/nc),(index)%(nc)].set_title(str(length)+'-mers')
            index+=1
    
    plt.suptitle(title,fontsize=25)
    plt.savefig(r.outputDir+'inlen.png')

def trading_iteration_figure2():
    #file_handler = open("simulation_result/FBSDE_expXdrift_linarYdrift/" + 
    #                "model_type_expXdrift_linearYdrift_A0.1_M_space10.0_M_time3.0_beta0.1_delta_space0.005_delta_time0.01_iter_number6_kappa3.0_m0.1_num_MC_path100000_sigma2.0_x_02.0",
    #                'rb')
    file_handler = open("simulation_result/FBSDE_expXdrift_linarYdrift/" +                  
                     "model_type_expXdrift_linearYdrift_A0.05_M_space15.0_M_time3.0_beta0.1_delta_space0.005_delta_time0.01_iter_number6_kappa3.0_m0.1_num_MC_path100000_sigma2.0_x_02.0",
                     'rb')
    tmp = pickle.load(file_handler)
    _, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 10))
    plt.suptitle("convergence of outputs from each iteration",
            y=0.95, size=18)
    iteration_figure(tmp, 
                 [0,1,2,3,-2,-1],
                 figure_ax=ax1,
                 x_array_limit=3,
                 xlabel="time",
                 ylabel="function value",
                 other=ax1.plot(np.arange(0,3,0.01),
                        trading_linear_theo(x_0=2,
                        m=0.1,
                        A=0.1,
                        kappa=3,
                        beta=0.1,
                        t_space=np.arange(0,3,0.01)), "--",
                        label="linearization result"))
    iteration_difference_figure(tmp, 
                            [0,1,2,3,-2],
                            figure_ax=ax2,
                            x_array_limit=3,
                            xlabel="time",
                            ylabel="difference from final output",
                            other=ax2.plot(np.arange(0,3,0.01),
                        np.zeros(len(np.arange(0,3,0.01))), "black"))
    
    plt.savefig("trading_iteration_figure.png")

def signmag_plot(a, b, z, ref):
    imdata1 = np.sign(ref)
    cmap1 = plt.cm.RdBu
    cmap1.set_bad('k', 1)

    imdata2 = np.log10(np.abs(ref))
    cmap2 = plt.cm.YlOrRd
    cmap2.set_bad('k', 1)

    fig, axarr = plt.subplots(ncols=2, figsize=(12, 6))
    axarr[0].pcolormesh(a, b, imdata1, cmap=cmap1, vmin=-1, vmax=1)
    im = axarr[1].pcolormesh(a, b, imdata2, cmap=cmap2,
                             vmin=np.percentile(imdata2,  5),
                             vmax=np.percentile(imdata2, 95))

    for ax in axarr:
        ax.set_xlim((np.min(a), np.max(a)))
        ax.set_ylim((np.min(b), np.max(b)))
        ax.set_xlabel("a")
        ax.set_ylabel("b")
        ax.set(adjustable='box-forced', aspect='equal')

    fig.subplots_adjust(right=0.8)
    cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.7])
    fig.colorbar(im, cax=cbar_ax)

    axarr[0].set_title("Sign of hyp1f1")
    axarr[1].set_title("Magnitude of hyp1f1")
    plt.suptitle("z = {:.2e}".format(np.float64(z)))

    return fig

def plot_forces_violinplots(experiment):
    ensemble = experiment.observations.kinematics
    ensembleF = ensemble.loc[
        (ensemble['position_x'] > 0.25) & (ensemble['position_x'] < 0.95),
        ['totalF_x', 'totalF_y', 'totalF_z',
         'randomF_x', 'randomF_y', 'randomF_z',
         'upwindF_x',
         'wallRepulsiveF_x', 'wallRepulsiveF_y', 'wallRepulsiveF_z',
         'stimF_x', 'stimF_y', 'stimF_z']] #== Nans
    # plot Forces
    #    f, axes = plt.subplots(2, 2, figsize=(9, 9), sharex=True, sharey=True)
    ##    forcefig = plt.figure(5, figsize=(9, 8))
    ##    gs2 = gridspec.GridSpec(2, 2)
    ##    Faxs = [fig.add_subplot(ss) for ss in gs2]
    forcefig = plt.figure()
    #    Faxs1 = forcefig.add_subplot(211)
    #    Faxs2 = forcefig.add_subplot(212)
    sns.violinplot(ensembleF, lw=3, alpha=0.7, palette="Set2")
    #    tF = sns.jointplot('totalF_x', 'totalF_y', ensemble, kind="hex", size=10)
    plt.suptitle("Force distributions")
    #    plt.xticks(range(4,((len(alignments.keys())+1)*4),4), [i[1] for i in medians_sgc], rotation=90, fontsize = 4)
    plt.tick_params(axis='x', pad=4)
    plt.xticks(rotation=40)
    #    remove_border()
    plt.tight_layout(pad=1.8)
    plt.ylabel("Force magnitude distribution (newtons)")

    fileappend, path, agent = get_agent_info(experiment.agent)

    plt.savefig(os.path.join(path, "Force Distributions" + fileappend + FIG_FORMAT))
    plt.show()