def plot_2D_offsets(df, stat_key, prefices=["Delta_x_", "Delta_y_"],
                    suffices=["_x", "_y"], prefix_bool=1):
    for i, key in enumerate(stat_key) :
        plt.axes().set_aspect('equal')
        if prefix_bool:
            fixed_keys = [prefices[0] + key, prefices[1] + key]
        else:
            fixed_keys = [key + suffices[0], key + suffices[1]]

        plt.plot(
            np.array(df[fixed_keys[0]]),
            np.array(df[fixed_keys[1]]), 'b.', alpha=0.05
            )
        biweight_loc = (
            biweight_location(df[fixed_keys[0]]),
            biweight_location(df[fixed_keys[1]]))

        # The red cross is the biweight location along each dimension
        plt.plot(biweight_loc[0], biweight_loc[1],
             'rx', mew=2.)
        plt.tick_params(labeltop='off', labelright='off')
        plt.axes().yaxis.set_ticks_position('left')
        plt.axes().xaxis.set_ticks_position('bottom')
        plt.xlim(-300, 300)
        plt.ylim(-300, 300)
        plt.title(key + ', biweight_loc = {0:.2f}, {1:.2f}'.format(
            *biweight_loc))

        plt.show()
        plt.clf()

def plot_images(data_list, data_shape="auto", fig_shape="auto"):
    """
    plotting data on current plt object.
    In default,data_shape and fig_shape are auto.
    It means considered the data as a sqare structure.
    """
    n_data = len(data_list)
    if data_shape == "auto":
        sqr = int(n_data ** 0.5)
        if sqr * sqr != n_data:
            data_shape = (sqr + 1, sqr + 1)
        else:
            data_shape = (sqr, sqr)
    plt.figure(figsize=data_shape)

    for i, data in enumerate(data_list):
        plt.subplot(data_shape[0], data_shape[1], i + 1)
        plt.gray()
        if fig_shape == "auto":
            fig_size = int(len(data) ** 0.5)
            if fig_size ** 2 != len(data):
                fig_shape = (fig_size + 1, fig_size + 1)
            else:
                fig_shape = (fig_size, fig_size)
        Z = data.reshape(fig_shape[0], fig_shape[1])
        plt.imshow(Z, interpolation="nearest")
        plt.tick_params(labelleft="off", labelbottom="off")
        plt.tick_params(axis="both", which="both", left="off", bottom="off", right="off", top="off")
        plt.subplots_adjust(hspace=0.05)
        plt.subplots_adjust(wspace=0.05)

 def __init__(self):
     # 図を描く大きさと、図の変数名を宣言
     fig = plt.figure(figsize=(5, 5))
     ax = plt.gca()
      
     # 赤い壁を描く
     plt.plot([1, 1], [0, 1], color='red', linewidth=2)
     plt.plot([1, 2], [2, 2], color='red', linewidth=2)
     plt.plot([2, 2], [2, 1], color='red', linewidth=2)
     plt.plot([2, 3], [1, 1], color='red', linewidth=2)
      
     # 状態を示す文字S0～S8を描く
     plt.text(0.5, 2.5, 'S0', size=14, ha='center')
     plt.text(1.5, 2.5, 'S1', size=14, ha='center')
     plt.text(2.5, 2.5, 'S2', size=14, ha='center')
     plt.text(0.5, 1.5, 'S3', size=14, ha='center')
     plt.text(1.5, 1.5, 'S4', size=14, ha='center')
     plt.text(2.5, 1.5, 'S5', size=14, ha='center')
     plt.text(0.5, 0.5, 'S6', size=14, ha='center')
     plt.text(1.5, 0.5, 'S7', size=14, ha='center')
     plt.text(2.5, 0.5, 'S8', size=14, ha='center')
     plt.text(0.5, 2.3, 'START', ha='center')
     plt.text(2.5, 0.3, 'GOAL', ha='center')
      
     # 描画範囲の設定と目盛りを消す設定
     ax.set_xlim(0, 3)
     ax.set_ylim(0, 3)
     plt.tick_params(axis='both', which='both', bottom=False, top=False,
                     labelbottom=False, right=False, left=False, labelleft=False)
      
     self.fig = fig
     self.ax = ax

def chartProperties(counter,path):

    seen_properties = sorted(counter, key=lambda x: x[1],reverse=True)
    seen_values_pct = map(itemgetter(1), tupleCounts2Percents(seen_properties))
    seen_values_pct = ['{:.1%}'.format(item)for item in seen_values_pct]

    plt.figure()

    numberchart = plt.bar(range(len(seen_properties)), map(itemgetter(1), seen_properties), width=0.9,alpha=0.6)
    plt.xticks(range(len(seen_properties)), map(itemgetter(0), seen_properties),rotation=90,ha='left')

    plt.ylabel('Occurrences')

    plot_margin = 1.15
    x0, x1, y0, y1 = plt.axis()
    plt.axis((x0,
              x1,
              y0,
              y1*plot_margin))

    plt.tick_params(axis='both', which='major', labelsize=8)
    plt.tick_params(axis='both', which='minor', labelsize=8)
    plt.tight_layout()

    autolabel(numberchart,seen_values_pct)

    plt.savefig(path)
    plt.clf()

def test_rotated_labels_parameters_no_ticks():
    fig, ax = __plot()

    ax.xaxis.set_ticks([])

    plt.tick_params(axis='x',
                    which='both',
                    bottom='off',
                    top='off')
    plt.tick_params(axis='y',
                    which='both',
                    left='off',
                    right='off')

    # convert to tikz file
    _, tmp_base = tempfile.mkstemp()
    tikz_file = tmp_base + '_tikz.tex'

    matplotlib2tikz.save(
        tikz_file,
        figurewidth='7.5cm'
        )

    # close figure
    plt.close(fig)

    # delete file
    os.unlink(tikz_file)
    return

def plot_first_k_numbers(X,k):
    m = X.shape[0]
    k = min(m,k)
    j = int(round(k / 10.0))

    fig, ax = plt.subplots(j,10)

    for i in range(k):
        w=X[i,:]

        w=w.reshape(28,28)
        ax[i/10, i%10].imshow(w, cmap=plt.cm.gist_yarg, interpolation='nearest', aspect='equal')
        ax[i/10, i%10].axis('off')

    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom='off',      # ticks along the bottom edge are off
        top='off',         # ticks along the top edge are off
        labelbottom='off')
    plt.tick_params(
        axis='y',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        left='off',
        right='off',    # ticks along the top edge are off
        labelleft='off')

    fig.show()

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()

def graph(csv_file, filename, bytes2str):
    '''Create a line graph from a two column csv file.'''

    unit = configs['unit']
    date, value = np.loadtxt(csv_file, delimiter=',', unpack=True,
                             converters={0: bytes2str}
                             )
    fig = plt.figure(figsize=(10, 3.5))
    fig.add_subplot(111, axisbg='white', frameon=False)
    rcParams.update({'font.size': 9})
    plt.plot_date(x=date, y=value, ls='solid', linewidth=2, color='#FB921D',
                  fmt=':'
                  )
    title = "Sump Pit Water Level {}".format(time.strftime('%Y-%m-%d %H:%M'))
    title_set = plt.title(title)
    title_set.set_y(1.09)
    plt.subplots_adjust(top=0.86)

    if unit == 'imperial':
        plt.ylabel('inches')
    if unit == 'metric':
        plt.ylabel('centimeters')

    plt.xlabel('Time of Day')
    plt.xticks(rotation=30)
    plt.grid(True, color='#ECE5DE', linestyle='solid')
    plt.tick_params(axis='x', bottom='off', top='off')
    plt.tick_params(axis='y', left='off', right='off')
    plt.savefig(filename, dpi=72)

def stationarity(timeseries):
	#Determing rolling statistics
	rol_mean = timeseries.rolling(window=12).mean()
	rol_std = timeseries.rolling(window=12).std()
	#Plot rolling statistics:
	fig, ax = plt.subplots()
	plt.grid(color='grey', which='major', axis='y', linestyle='--')
	plt.plot(timeseries, color='blue', label='Original', linewidth=1.25)
	plt.plot(rol_mean, color='red', label='Rolling Mean', linewidth=1.25)
	plt.plot(rol_std, color='black', label = 'Rolling Std', linewidth=1.25)
	plt.legend(loc='best')
	title = headers[1], data[index].iloc[0], '-' ,data[index].iloc[-1]
	plt.title(title)
	plt.tick_params(axis="both", which="both", bottom="on", top="off",    
		            labelbottom="on", left="off", right="off", labelleft="on")
	ax.spines['right'].set_visible(False)
	ax.spines['top'].set_visible(False)
	ax.xaxis.set_ticks_position('bottom') 
	fig.title = ('stationarity.png')
	fig.savefig(fig.title, bbox_inches="tight")
	#Perform Dickey-Fuller test:
	print ('Results of Dickey-Fuller Test:\n')
	df_test = adfuller(timeseries, autolag='AIC')
	df_output = pd.Series(df_test[0:4], index=['Test Statistic','p-value','#Lags Used','No. of Observations Used'])
	for key,value in df_test[4].items():
	    df_output['Critical Value (%s)'%key] = value
	print (df_output.round(3))

def SetAxes(legend=False):
    f_b = 0.164
    f_star = 0.01
    err_b = 0.006
    err_star = 0.004
    f_gas = f_b - f_star
    err_gas = np.sqrt(err_b**2 + err_star**2)

    plt.axhline(y=f_gas, ls='--', c='k', label='', zorder=-1)
    x = np.linspace(.0,2.,1000)
    plt.fill_between(x, y1=f_gas - err_gas, y2=f_gas + err_gas, color='k', alpha=0.3, zorder=-1)
    plt.text(.6, f_gas+0.006, r'f$_{gas}$', verticalalignment='bottom', size='large')
    plt.xlabel(r'r/r$_{vir}$', size='x-large')
    plt.ylabel(r'f$_{gas}$ ($<$ r)', size='x-large')

    plt.xscale('log')
    plt.xticks([1./1.9, 1.33/1.9, 1, 1.5, 2.],[r'r$_{500}$', r'r$_{200}$', 1, 1.5, 2], size='large')
    #plt.yticks([.1, .2], ['0.10', '0.20'])
    plt.tick_params(length=10, which='major')
    plt.tick_params(length=5, which='minor')
    plt.xlim([0.4,1.5])
    plt.minorticks_on()

    if legend:
        plt.legend(loc=0, prop={'size':'small'}, markerscale=0.7, numpoints=1, ncol=2)

def histogram(x, y, xlabel=None, ylabel=None, title=None, out=None, highlight=None):
    plt.bar(range(len(x)), y, color='b', alpha=0.6, linewidth=0, align='center')

    # if highlight is not None:
    #     barlist[highlight].set_color('r')
    #     rect = barlist[highlight]
    #     height = rect.get_height()
    #     width = rect.get_width()
    #     plt.text(rect.get_x() + width/2., height+1, "%.2f%%" % float(height), ha='center', va='center')

    plt.xticks(range(len(x)), x, ha='center', va='top', rotation='vertical')
    plt.tick_params(axis='y', labelsize='x-small')
    # plt.xlim([-1, len(x)])
    # plt.ylim([0, y[0]+((y[0]/100)*10)])
    # plt.axis('auto')
    # plt.title(title)
    plt.gca().yaxis.grid(True)


    if xlabel != None and ylabel != None:
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)

    if out == None:
        plt.show()
    else:
        plt.savefig("../PLOT/" + out + ".jpg")
        plt.show()

def drawMean(ctype, filename, titulo, yaxis, comp, *args):
    plt.switch_backend('Qt4Agg')
    fig = plt.figure(figsize=(28, 5))
    ax = fig.add_subplot(1,2,1)
    means = []
    for v in comp:
        means.append(np.mean(v))
    col = ['r', 'b', 'g', 'm', 'c', 'y', 'r', 'b', 'g', 'm', 'c', 'y']
    if len(comp) == 16:
        col = ['r', 'r', 'g', 'g', 'g', 'b', 'b', 'b', 'c', 'c', 'c', 'm', 'm', 'm', 'y', 'y', 'y' ]
    if len(comp) == 12:
        col = ['r', 'r' ,'b','b','g', 'g','m', 'm', 'c', 'c', 'y' ,'y']
    if len(comp) == 13:
        col = ['r', 'r', 'r','b','b','g', 'g','m', 'm', 'c', 'c', 'y' ,'y']
    if len(comp) == 11:
        col = ['r' ,'b','b','g', 'g','m', 'm', 'c', 'c', 'y' ,'y']
    ax.grid( b=True, linestyle='-',  axis = 'y', linewidth=1, zorder=1) 
    if len(comp) < 7:
        ax.bar( range(1, int(ctype) + 1),  means, align='center', color=col[0: int(ctype) + 1], zorder = 10, width = 0.6)
    else:
        ax.bar( range(1, int(ctype) + 1),  means, align='center', color=col[0: int(ctype) + 1], zorder = 10)
    plt.title(titulo)
    labels = [i for i in args]
    plt.ylabel(yaxis+' (promedio)')
    plt.tick_params(labelsize = 7)
    plt.xticks(range(1, (int(ctype)) + 1), labels)
    plt.savefig(filename+' (promedio).png', bbox_inches = 'tight')
    plt.xticks()
    val_strings = ["%10.10f" % x for x in means ]
    with open(filename + '_promedio.txt', 'w') as file_:
        for k in range(0, len(comp)):
            file_.write(labels[k] + " " + val_strings[k] + "\n")

def histogram(X, z, colors=['#1F77B4', '#FF7F0E'], fname='plot.pdf',
              xlim=None, bins=500):
    """Plot histograms of 1-dimensional data.
    
    Input: X = data matrix
           z = class labels
           color = color for each class
           fname = output file

    Output: None

    """
    z_unique = np.unique(z)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    colors = iter(colors)
    for k in z_unique:
        idx = np.where(z==k)[0]
        x = X[idx]
        ax.hist(x, bins=bins, facecolor=next(colors), 
                histtype='stepfilled',
                alpha=.8, normed=1, linewidth=0.5)
    ax.set_xlabel(r'$x$')
    if xlim:
        ax.set_xlim(xlim)
    plt.tick_params(top='off', bottom='on', left='off', right='off',
            labelleft='off', labelbottom='on')
    for i, spine in enumerate(plt.gca().spines.values()):
        if i !=2:
            spine.set_visible(False)
    frame = plt.gca()
    frame.axes.get_yaxis().set_visible(False)
    fig.tight_layout()
    fig.savefig(fname)

def draw_filters(W, cols=20, fig_size=(10, 10), filter_shape=(28, 28),
                 filter_standardization=False):
    border = 2
    num_filters = len(W)
    rows = int(np.ceil(float(num_filters) / cols))
    filter_height, filter_width = filter_shape

    if filter_standardization:
        W = preprocessing.scale(W, axis=1)
    image_shape = (rows * filter_height + (border * rows),
                   cols * filter_width + (border * cols))
    low, high = W.min(), W.max()
    low = (3 * low + high) / 4
    high = (low + 3 * high) / 4
    all_filter_image = np.random.uniform(low=low, high=high,
                                         size=image_shape)
    all_filter_image = np.full(image_shape, W.min(), dtype=np.float32)

    for i, w in enumerate(W):
        start_row = (filter_height * (i / cols) +
                     (i / cols + 1) * border)
        end_row = start_row + filter_height
        start_col = (filter_width * (i % cols) +
                     (i % cols + 1) * border)
        end_col = start_col + filter_width
        all_filter_image[start_row:end_row, start_col:end_col] = \
            w.reshape(filter_shape)

    plt.figure(figsize=fig_size)
    plt.imshow(all_filter_image, cmap=plt.cm.gray,
               interpolation='none')
    plt.tick_params(axis='both',  labelbottom='off',  labelleft='off')
    plt.show()

def plot_nice(data, xlabel = 'Time (s)',ylabel='ylabel',title='Title',color='b', dims = 1, filename='nice',
              size=(10,5), fontsize= 10,lw=2,marker='',xticks=[],yticks=[]):
    fig = plt.figure()
    ax = plt.gca()
    fig.set_size_inches(size)
    if dims == 1:
        plt.plot(data,linewidth = lw, color =color ,marker=marker)
    else:
        plt.plot(data[0],data[1], linewidth = lw, color=color)
    if not ARTISTIC_PLOT:
        plt.xlabel(xlabel,fontsize=fontsize)
    plt.ylabel(ylabel,fontsize=fontsize)
    if SHOW_TITLE:
        plt.title(title,fontsize=fontsize)
    plt.tick_params(labelsize=fontsize)
    ax.spines['top'].set_color('none')
    ax.spines['right'].set_color('none')
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    plt.gcf().subplots_adjust(bottom=.4,left=.2)

    if ARTISTIC_PLOT:
        ax.spines['bottom'].set_color('none')
    if len (xticks ) > 0 and not ARTISTIC_PLOT:
        plt.xticks(xticks,xlabels)
    #if len (yticks ) > 0:
    plt.yticks(yticks)
    if ARTISTIC_PLOT:
        ax.spines['bottom'].set_color('none')
        plt.xticks([])
    if SAVE_PLOTS:
        plt.savefig('results/' + filename+'.eps')
        plt.savefig('results/' + filename+'.png')

def ARMA(timeseries):
	# Input how many lags
	print ('\nHow many AR lags?')
	p = int(input())
	print ('\nHow many MA lags?')
	q = int(input())
	print ('\nDifferenced? Y/N')
	d = 1 if input() == 'Y' else 0
	arma_mod = sm.tsa.ARIMA(timeseries, order=(p, d, q))
	arma_res = arma_mod.fit(trend='nc', disp=-1)
	print (arma_res.summary())
	# Plot ARIMA w/ predictions
	fig, ax = plt.subplots()
	ax = timeseries.ix[1:].plot(ax=ax)
	plt.tick_params(axis="both", which="both", bottom="on", top="off",    
		                labelbottom="on", left="on", right="off", labelleft="on")
	print('When would you like to start forecast? DD-MM-YYY')
	pred_start = (input())
	print('When would you like to stop forecast? DD-MM-YYY (e.g. one year ahead)')
	pred_end = (input())
	arma_res.plot_predict(start=pred_start, end=pred_end, dynamic=False, ax=ax, plot_insample=False)	
	plt.xlabel("Year", fontsize=14) 
	plt.title('ARIMA', fontsize=15)
	plt.tight_layout()
	fig.savefig('arma.png', bbox_inches="tight")

def plot_2hist(false_values, true_values, hist_bins=None, labels=None, hist_title=None ):
    
    if not hist_bins:
        try:
            hist_bins = ( min( min(false_values), min(true_values)), max( max(false_values), max(true_values)), 50 )
        except ValueError:
            hist_bins = (0, 1, 2)
        
    bins = numpy.linspace( hist_bins[0], hist_bins[1], hist_bins[2] )
    
    if not labels:
        label0 = namestr(false_values)
        label1 = namestr(true_values)
    else:
        label0 = labels[0]
        label1 = labels[1]
    
    pyplot.hist([ [false_values], [true_values] ], bins, color=('red', 'blue'), label=(label0,label1), normed=True )
    
    pyplot.legend(loc='upper right', fontsize=18)
    pyplot.title( hist_title,  fontsize=24)
    
    pyplot.tick_params(axis='x', labelsize=16)
    pyplot.tick_params(axis='y', labelsize=16)

    return 1

def plot_date_bars(bin_data, bin_edges, title, ylabel, fname):
    """
    Semi-generic function to plot a bar graph, x-label is fixed to "date" and the
    x-ticks are formatted accordingly.

    To plot a histogram, the histogram data must be calculated manually outside
    this function, either manually or using :py:func`numpy.histogram`.

    :param bin_data: list of data for each bin
    :param bin_edges: list of bin edges (:py:class:`datetime.date` objects), its
                      length must be ``len(data)+1``
    :param title: title of the plot
    :param ylabel: label of y-axis
    :param fname: output file name
    """
    import matplotlib.pyplot as plt
    from matplotlib.dates import date2num, num2date
    from matplotlib import ticker

    plt.figure()  # clear previous figure
    plt.title(title)
    plt.xlabel("date")
    plt.ylabel(ylabel)

    # plot the bars, width of the bins is assumed to be fixed
    plt.bar(date2num(bin_edges[:-1]), bin_data, width=date2num(bin_edges[1]) - date2num(bin_edges[0]))

    # x-ticks formatting
    plt.gca().xaxis.set_major_formatter(ticker.FuncFormatter(lambda numdate, _: num2date(numdate).strftime('%Y-%m-%d')))
    plt.gcf().autofmt_xdate()
    plt.tick_params(axis="x", which="both", direction="out")
    plt.xticks([date2num(ts) for ts in bin_edges if ts.month % 12 == 1])

    plt.savefig(fname, papertype="a4")

def makeSigBkg(all_outputs, targets, label,
              dir='/afs/cern.ch/user/j/jpavezse/systematics',model_g='mlp',
              print_pdf=False,legends=None, title=''):
  '''
  make plots for ROC curve of classifier and
  test data.
  '''
  tprs = []
  fnrs = []
  aucs = []
  thresholds = np.linspace(0,1.0,150) 
  fig = plt.figure()
  for k,(outputs,target) in enumerate(zip(all_outputs,targets)):
    fnrs.append(np.array([float(np.sum((outputs > tr) * (target == 0)))/float(np.sum(target == 0)) for tr in thresholds]))
    fnrs[-1] = fnrs[-1].ravel()
    tprs.append(np.array([float(np.sum((outputs < tr) * (target == 1)))/float(np.sum(target == 1)) for tr in thresholds]))
    tprs[-1] = tprs[-1].ravel()
    aucs.append(auc(tprs[-1],fnrs[-1]))
    plt.plot(tprs[-1], fnrs[-1], label='ROC {0} (area = {1:.2f})'.format(
       legends[k],aucs[-1]))
  plt.xlim([0.0, 1.0])
  plt.ylim([0.0, 1.05])
  plt.xlabel('Signal Efficiency',fontsize=11)
  plt.ylabel('Background Rejection',fontsize=11)
  plt.tick_params(axis='both', labelsize=10)
  plt.title('{0}'.format(title))
  plt.legend(loc="lower left",frameon=False, fontsize=11)
  #np.savetxt('{0}/plots/{1}/results/{2}.txt'.format(dir,model_g,label),np.column_stack((tpr,fnr)))
  plt.savefig('{0}/plots/{1}/{2}.png'.format(dir,model_g,label))
  if print_pdf == True:
    plt.savefig('{0}/plots/{1}/{2}.pdf'.format(dir,model_g,label))
  plt.close(fig)
  plt.clf()

 def plot_sleipner_thick_contact(self, years, gwc = False, sim_title = ''):
     if gwc == True:
         tc_str = 'contact'
     else:
         tc_str = 'thickness'
     yr_indices = self.get_plan_year_indices(years)
     size = 14
     font = {'size' : size}
     matplotlib.rc('font', **font)
     fig = plt.figure(figsize=(10.0, 2.5), dpi = 960)
     middle = len(years) * 10
     pos = 100 + middle
     for n in range(len(yr_indices)):
         pos +=1
         ax = fig.add_subplot(pos)
         xf = []
         yf = []
         kf = []
         for i in range(self.nx):
             tempx = []
             tempy = []
             tempk = []
             for j in range(self.ny):
                 x = self.x[(i, j, 0)]
                 y = self.y[(i, j, 0)]
                 tn = yr_indices[n]
                 thick, contact = self.get_thick_contact(i, j, tn)
                 tempx.append(x)
                 tempy.append(y)
                 if gwc == True:
                     tempk.append(contact)
                 else:
                     tempk.append(thick)
             xf.append(tempx)
             yf.append(tempy)
             kf.append(tempk)
         xp = np.asarray(xf)
         yp = np.asarray(yf)
         kp = np.asarray(kf)
         N = 10
         contour_label = False
         ax_label = False
         c = ax.contourf(xp, yp, kp, N)
         plt.tick_params(which='major', length=3, color = 'w')
         if n == len(years) - 1:
             fig.subplots_adjust(right=0.84)
             cb_axes = fig.add_axes([0.85, 0.15, 0.05, 0.7])
             plt.tick_params(which='major', length=3, color = 'k')
             cb = fig.colorbar(c, cax = cb_axes, format = '%.2f')
             cb.set_ticks(np.linspace(np.amin(kp), np.amax(kp), N))
             cb.set_label(tc_str + ': [m]')
         if n != 0:
             ax.set_yticklabels([])
         ax.set_xticklabels([])
         ax.set_title(str(years[n]))
         ax.axis([0, 3000, 0, 6000])
         ax.xaxis.set_ticks(np.arange(0,3500,1000))
     plt.savefig(sim_title + '_' +  tc_str + '.pdf', fmt = 'pdf')
     plt.clf()
     return 0