Python matplotlib.pyplot.hlines() Examples

def plot(args, timings):
    for name, cls_timings in timings:
        xs, relative_timings = zip(*cls_timings)
        ys = [r.normal / r.compiled for r in relative_timings]
        for x, y in zip(xs, ys):
            print(xs, ys)
        if args.show_plot:
            plt.plot(xs, ys, '-o', label=name)
            plt.hlines(1.0, np.min(xs), np.max(xs), 'k')

    if not args.show_plot:
        return

    plt.xlabel('Number of weak learners')
    plt.ylabel('Relative speedup')
    plt.axis('tight')
    plt.legend()
    plt.gca().set_ylim(bottom=0)
    title, suptitle = titles(args)
    plt.title(title)
    plt.suptitle(suptitle, fontsize=3)
    filename = "timings{0}.png".format(hash(str(args)))
    plt.savefig(filename, dpi=72) 

def flip_plot(min_exp, max_exp):
    """
    Assumes min_exp and min_exp positive integers; min_exp < max_exp
    Plots results of 2**min_exp to 2**max_exp coin flips
    抛硬币的次数为2的min_exp次方到2的max_exp次方
    一共进行了 2**max_exp - 2**min_exp 轮实验，每轮实验抛硬币次数逐渐增加
    """

    ratios = []
    x_axis = []
    for exp in range(min_exp, max_exp + 1):
        x_axis.append(2**exp)
    for numFlips in x_axis:
        num_heads = 0  # 初始化，硬币正面朝上的计数为0
        for n in range(numFlips):
            if random.random() < 0.5:  # random.random()从[0, 1)随机的取出一个数
                num_heads += 1  # 当随机取出的数小于0.5时，正面朝上的计数加1
        num_tails = numFlips - num_heads  # 得到本次试验中反面朝上的次数
        ratios.append(num_heads/float(num_tails))  # 正反面计数的比值
    plt.title('Heads/Tails Ratios')
    plt.xlabel('Number of Flips')
    plt.ylabel('Heads/Tails')
    plt.plot(x_axis, ratios)
    plt.hlines(1, 0, x_axis[-1], linestyles='dashed', colors='r')
    plt.show() 

def view_palette(*args):
    if len(args) > 1:
        f, ax = plt.subplots(1, len(args), figsize=(3 * len(args), 3))
        for i, name in enumerate(args):
            check_key(name)
            cycle = palettes[name]
            for j, c in enumerate(cycle):
                ax[i].hlines(j, 0, 1, colors=c, linewidth=15)
            ax[i].set_title(name)
            despine(ax[i], True)
        plt.show()
    elif len(args) == 1:
        f = plt.figure(figsize=(3, 3))
        check_key(args[0])
        cycle = palettes[args[0]]
        for j, c in enumerate(cycle):
            plt.hlines(j, 0, 1, colors=c, linewidth=15)
        plt.title(args[0])
        despine(plt.axes(), True)
        f.tight_layout()
        plt.show()
    else:
        raise NotImplementedError("ERROR: supply a palette to plot. check vapeplot.available() for available palettes") 

def rfReg(X_rfReg,y_rfReg,predFeat=False):     
     X=X_rfReg
     y=y_rfReg
     X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.4,random_state=1)
     forest=RandomForestRegressor(n_estimators=100,criterion='mse',random_state=1,n_jobs=-1)
     forest.fit(X_train,y_train)
     y_train_pred=forest.predict(X_train)
     y_test_pred=forest.predict(X_test)
     print('MSE train:%.3f,test:%.3f'%(mean_squared_error(y_train,y_train_pred),mean_squared_error(y_test,y_test_pred)))
     print('R^2 train:%.3f,test:%.3f'%(r2_score(y_train,y_train_pred),r2_score(y_test,y_test_pred)))
    
     plt.scatter(y_train_pred,y_train_pred-y_train,c='blue',marker='o',label='Training data')
     plt.scatter(y_test_pred,y_test_pred-y_test,c='lightgreen',marker='s',label='Test data')
     plt.xlabel('Predicted values')
     plt.ylabel('Residuals')
     plt.legend(loc='upper left')
     plt.hlines(y=0,xmin=0,xmax=6,lw=2,color='red')
     plt.xlim([0,6])
     plt.show()
     
     if type(predFeat).__module__=='numpy':
         return forest.predict(predFeat) 

def LR_multiRegression(X_lrm,y_lrm,predFeat=False):    
    X=X_lrm
    y=y_lrm
    X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.3,random_state=0)
    slr=LinearRegression()
    slr.fit(X_train,y_train)
    y_train_pred=slr.predict(X_train)
    y_test_pred=slr.predict(X_test)
    
    plt.scatter(y_train_pred,y_train_pred-y_train,c='blue',marker='o',label='Training data')
    plt.scatter(y_test_pred,y_test_pred-y_test,c='lightgreen',marker='s',label='Test data')
    plt.xlabel('Predicted values')
    plt.ylabel('Residuals')
    plt.legend(loc='upper left')
    plt.hlines(y=0,xmin=-10,xmax=7,lw=2,color='red')
    plt.xlim([0,7])
    plt.show()
      
    print(slr.coef_,slr.intercept_)
    print('MSE train:%.3f,test:%.3f'%(mean_squared_error(y_train,y_train_pred),mean_squared_error(y_test,y_test_pred))) #MSE， float or ndarray of floats,A non-negative floating point value (the best value is 0.0), or an array of floating point values, one for each individual target.
    print('R^2 train:%.3f,test:%.3f'%(r2_score(y_train,y_train_pred),r2_score(y_test,y_test_pred)))

    if type(predFeat).__module__=='numpy':
        return slr.predict(predFeat) 

def visualize_anomaly(y_true, reconstruction_error, threshold):
    error_df = pd.DataFrame({'reconstruction_error': reconstruction_error,
                             'true_class': y_true})
    print(error_df.describe())

    groups = error_df.groupby('true_class')
    fig, ax = plt.subplots()

    for name, group in groups:
        ax.plot(group.index, group.reconstruction_error, marker='o', ms=3.5, linestyle='',
                label="Fraud" if name == 1 else "Normal")

    ax.hlines(threshold, ax.get_xlim()[0], ax.get_xlim()[1], colors="r", zorder=100, label='Threshold')
    ax.legend()
    plt.title("Reconstruction error for different classes")
    plt.ylabel("Reconstruction error")
    plt.xlabel("Data point index")
    plt.show() 

def statistics(df, keys=None, plot=False):
    data = {}
    for key in keys:
        try:
            df_meas = new_sample(df[key])
            dt = find_meas_period(df_meas)
            stddev = np.sqrt(df_meas.var())
            mean = df_meas.mean()
            variance = stddev ** 2
            noise_power = variance * dt
        except Exception:
            dt = 0
            stddev = 0
            mean = 0
            noise_power = 0
            variance = 0
        data[key + '_stddev'] = stddev
        data[key + '_variance'] = variance
        data[key + '_mean'] = mean
        data[key + '_dt'] = dt
        data[key + '_noise_power'] = noise_power
        if plot:
            import matplotlib.pyplot as plt
            plt.figure()
            plt.title(key + ' statistics')
            df[key].plot(alpha=0.5)
            plt.hlines(mean, df.index[0], df.index[-1], 'k')
            plt.hlines([mean + stddev, mean - stddev],
                       df.index[0], df.index[-1], 'r')
    return data 

def edit_score(input_df, pred_df, filepath="edit_score.csv", plot=True):
    assert np.all(input_df['id'].values == pred_df['id'].values)
    if not np.all(np.array([len(x) for x in pred_df['expected']]) == np.array([len(x) for x in input_df['input']])):
        print("Some lengths do not match!")
        return None, None 
    output_df = input_df.copy().reset_index(drop=True)
    lev_dist = [Levenshtein.distance(x, y) for x, y in zip(input_df['expected'], pred_df['expected'])]
    ham_dist = [ham_distance(x, y) for x, y in zip(input_df['expected'], pred_df['expected'])]
    lev_score = np.mean(lev_dist)
    ham_score = np.mean(ham_dist)

    total_ham = np.sum(ham_dist)
    total_len = input_df['expected'].map(len).sum()
    accuracy = 1 - total_ham / total_len

    output_df['predicted'] = pred_df['expected'].values
    output_df['levdist'] = np.array(lev_dist)
    output_df['hamdist'] = np.array(ham_dist)
    output_df['levpercent'] = output_df['levdist'] / output_df['len']
    output_df['hampercent'] = output_df['hamdist'] / output_df['len']
    output_df['accuracy'] = 1 - output_df['hampercent']
    ham_percent = np.mean(output_df['hampercent'])
    mean_acc = np.mean(output_df['accuracy'])

    output_df.to_csv(filepath, index=False)
    print_str = "total acc: {:.4f}, mean acc: {:.4f}, lev: {:.1f}, ham: {:.1f}".format(accuracy, mean_acc, lev_score, ham_score)
    print(print_str)
    output_df.plot("len", "accuracy", kind="scatter")
    plt.hlines(y=accuracy, xmin=0, xmax=output_df['len'].max())
    plt.title(print_str)
    plt.savefig(filepath.replace(".csv", "_plot.png"))
    if plot:
        plt.show()
    plt.close()
    return accuracy, output_df

# Computes and returns the n-grams of a particualr sequence, defaults to trigrams 

def plot_grad_flow(named_parameters,n_iter):

    '''Plots the gradients flowing through different layers in the net during training.
    Can be used for checking for possible gradient vanishing / exploding problems.

    Usage: Plug this function in Trainer class after loss.backwards() as
    "plot_grad_flow(self.model.named_parameters())" to visualize the gradient flow'''
    ave_grads = []
    max_grads = []
    layers = []
    for n, p in named_parameters:
        if (p.requires_grad) and ("bias" not in n):
            layers.append(n)
            ave_grads.append(p.grad.abs().mean())
            max_grads.append(p.grad.abs().max())
    plt.bar(np.arange(len(max_grads)), max_grads, alpha=0.1, lw=1, color="c")
    plt.bar(np.arange(len(max_grads)), ave_grads, alpha=0.1, lw=1, color="b")
    plt.hlines(0, 0, len(ave_grads) + 1, lw=2, color="k")
    plt.xticks(range(0, len(ave_grads), 1), layers, rotation="vertical")
    plt.xlim(left=0, right=len(ave_grads))
    plt.ylim(bottom=-0.001, top=0.02)  # zoom in on the lower gradient regions
    plt.xlabel("Layers")
    plt.ylabel("average gradient")
    plt.title("Gradient flow")
    plt.grid(True)
    plt.legend([Line2D([0], [0], color="c", lw=4),
                Line2D([0], [0], color="b", lw=4),
                Line2D([0], [0], color="k", lw=4)], ['max-gradient', 'mean-gradient', 'zero-gradient'])
    #plt.savefig('./model/pred/Grad_Flow_' + str(n_iter - 1)) 

def simpleLR(X0,X1,w0,w1,w2):
    y=w0+w1*X0  #单解释变量
    plt.scatter(X0,y,c='r',marker='x')
    plt.plot(X0,y,label='linear fit',linestyle='--')
    plt.hlines(y=0,xmin=-10,xmax=60,lw=1,color='red')
    plt.vlines(x=0,ymin=-100,ymax=600,lw=1,color='red')    
#    plt.plot(,linestyle='-')
    plt.show()
    
    if X1.shape and w1 and w2:
        y=w0+w1*X0+w2*X1  #多解释变量
        fig=plt.figure()
        ax=Axes3D(fig)
        ax.scatter(X0, X1, y)
        plt.show() 

def plot_baseline(algo, task, color='r'):
  """Plot baseline as a horizontal dashed line."""
  key = '%s.%s' % (algo, task)
  if key in BASELINES:
    baseline = BASELINES[key]
    xlims = plt.xlim()
    plt.hlines(baseline, xlims[0], xlims[1], linestyles='dashed', colors=color) 

def plot_lcurve(self, ax=None, guides=True):
        """
        Make a plot of the data-misfit x regularization values.

        The estimated corner value is shown as a blue triangle.

        Parameters:

        * ax : matplotlib Axes
            If not ``None``, will plot the curve on this Axes instance.
        * guides : True or False
            Plot vertical and horizontal lines across the corner value.

        """
        if ax is None:
            ax = mpl.gca()
        else:
            mpl.sca(ax)
        x, y = self.dnorm, self.mnorm
        if self.loglog:
            mpl.loglog(x, y, '.-k')
        else:
            mpl.plot(x, y, '.-k')
        if guides:
            vmin, vmax = ax.get_ybound()
            mpl.vlines(x[self.corner_], vmin, vmax)
            vmin, vmax = ax.get_xbound()
            mpl.hlines(y[self.corner_], vmin, vmax)
        mpl.plot(x[self.corner_], y[self.corner_], '^b', markersize=10)
        mpl.xlabel('Data misfit(data norm)')
        mpl.ylabel('Regularization(model norm)') 

def plot_multi_pr_curves(plot_tuples, plot_title='Precision Recall Curves',
                         figsize=(12, 8), xlim=(0, 1), ylim=(0, 1),
                         basic_font_size=14):
    plt.figure(figsize=figsize)

    for eval_infos, line_name, line_color in plot_tuples:
        precs = eval_infos[0]
        recalls = eval_infos[1]
        avg_prec = eval_infos[3]
        f1_score = eval_infos[6]
        plt.step(recalls, precs,
                 label=line_name + ' (AUC {0:.3f}, F1 {1:.3f})'.format(avg_prec, f1_score),
                 color=line_color)

        dec_prec = eval_infos[4]
        dec_recall = eval_infos[5]
        plt.plot(dec_recall, dec_prec, 'o', color=line_color, markersize=8)
        plt.vlines(dec_recall, 0, dec_prec, linestyles='dashed', colors=line_color)
        plt.hlines(dec_prec, 0, dec_recall, linestyles='dashed', colors=line_color)

    plt.legend(fontsize=basic_font_size)
    plt.title(plot_title, fontsize=basic_font_size+ 2)
    plt.xlabel('Recall', fontsize=basic_font_size)
    plt.ylabel('Precision', fontsize=basic_font_size)
    plt.xticks(fontsize=basic_font_size)
    plt.yticks(fontsize=basic_font_size)
    plt.xlim(xlim)
    plt.ylim(ylim) 

def plot_dendrogram(a_timer: PipelineTimer,
                    linkage_matrix,
                    threshold: float,
                    vehicle_number: str,
                    force: bool = False):
    dendrogram_filename = "dendrogram_" + vehicle_number + "." + figure_format
    if path.isfile(dendrogram_filename):
        if force:
            remove(dendrogram_filename)
        else:
            print("Dendrogram already plotted. Skipping...")
            return
    plt.figure(figsize=(7, 7), dpi=600)
    R: dict = dendrogram(Z=linkage_matrix, orientation='top', distance_sort='ascending', no_labels=True)
    plt.title("Dendrogram of Agglomerative Clustering for Vehicle " + vehicle_number)
    plt.xlabel("Signals Observed")
    plt.ylabel("Single Linkage Cluster Merge Distance")
    xmin, xmax = plt.xlim()
    # Add a 25% opacity dashed black line to the entropy gradient plot at one boundary of each sub-flow
    plt.hlines(y=threshold, xmin=xmin, xmax=xmax, alpha=0.25, colors='black', linestyle='dashed',
               label='cluster threshold')
    plt.legend(loc='upper right')

    print("\tPlotting dendrogram and saving to " + dendrogram_filename)

    savefig(dendrogram_filename,
            bbox_iches='tight',
            pad_inches=0.0,
            dpi=600,
            format=figure_format,
            transparent=figure_transp)
    plt.close()
    print("\t\tComplete...") 

def add_plot_classification(fig, df, depth_max, depth_min, title, num_col, z_NAP):
    """
    Add to the plot the selected classification.

    :param fig: Original figure.
    :param depth_max: Maximum depth.
    :param depth_min: Minimum depth.
    :return:
    """
    ax = fig.add_subplot(1, num_col, 3)
    df = df.copy()
    df["soil_type"].loc[df["soil_type"].isna()] = "UNKNOWN"
    for st in np.unique(df["soil_type"]):
        partial_df = df[df["soil_type"] == st]
        if z_NAP:
            plt.hlines(
                y=partial_df["elevation_with_respect_to_NAP"],
                xmin=0,
                xmax=1,
                colors=partial_df["colour"],
                label=st,
            )
        else:
            plt.hlines(
                y=partial_df["depth"],
                xmin=0,
                xmax=1,
                colors=partial_df["colour"],
                label=st,
            )
    ax.set_xticks([])
    ax.set_title(f"{title}")
    plt.ylim(depth_max, depth_min)
    return fig 

def test_tributary(self, default_calving):

        _, ds1, df_diag1 = default_calving

        model = FluxBasedModel(bu_tidewater_bed(split_flowline_before_water=5),
                               mb_model=ScalarMassBalance(),
                               is_tidewater=True, calving_use_limiter=True,
                               smooth_trib_influx=False,
                               flux_gate=[0.06, 0], do_kcalving=True,
                               calving_k=0.2)
        _, ds2 = model.run_until_and_store(3000)
        df_diag2_a = model.get_diagnostics(fl_id=0)
        df_diag2_b = model.get_diagnostics(fl_id=1)

        assert_allclose(model.volume_m3 + model.calving_m3_since_y0,
                        model.flux_gate_m3_since_y0)
        assert_allclose(ds2.calving_m3[-1], model.calving_m3_since_y0)
        assert_allclose(ds2.volume_bsl_m3[-1], model.volume_bsl_km3 * 1e9)
        assert_allclose(ds2.volume_bwl_m3[-1], model.volume_bwl_km3 * 1e9)

        # should be veeery close
        rtol = 5e-4
        assert_allclose(ds1.volume_m3[-1], ds2.volume_m3[-1], rtol=rtol)
        assert_allclose(ds1.calving_m3[-1], ds2.calving_m3[-1], rtol=rtol)
        assert_allclose(ds1.volume_bsl_m3[-1], ds2.volume_bsl_m3[-1],
                        rtol=rtol)
        assert_allclose(ds2.volume_bsl_m3, ds2.volume_bwl_m3, rtol=rtol)

        df_diag1['surface_h_trib'] = np.append(df_diag2_a['surface_h'],
                                               df_diag2_b['surface_h'])

        if do_plot:
            f, ax = plt.subplots(1, 1, figsize=(12, 5))
            df_diag1[['surface_h', 'surface_h_trib',
                      'bed_h']].plot(ax=ax, color=['C3', 'C1', 'k'])
            plt.hlines(0, 0, 60000, color='C0', linestyles=':')
            plt.ylim(-350, 800)
            plt.ylabel('Altitude [m]')
            plt.show() 

def visualize_reconstruction_error(reconstruction_error, threshold):
    plt.plot(reconstruction_error, marker='o', ms=3.5, linestyle='',
             label='Point')

    plt.hlines(threshold, xmin=0, xmax=len(reconstruction_error)-1, colors="r", zorder=100, label='Threshold')
    plt.legend()
    plt.title("Reconstruction error")
    plt.ylabel("Reconstruction error")
    plt.xlabel("Data point index")
    plt.show() 

def graph_barcode(data, ph, homology_group=0):
    persistence = ph.transform(data)
    # this function just produces the barcode graph for each homology group
    xstart = [s[1][0] for s in persistence if s[0] == homology_group]
    xstop = [s[1][1] for s in persistence if s[0] == homology_group]
    y = [0.1 * x + 0.1 for x in range(len(xstart))]
    plt.hlines(y, xstart, xstop, color='b', lw=4)
    # Setup the plot
    ax = plt.gca()
    plt.ylim(0, max(y) + 0.1)
    ax.yaxis.set_major_formatter(plt.NullFormatter())
    plt.xlabel('epsilon')
    plt.ylabel("Betti dim %s" % (homology_group,))
    plt.show() 

def plot_date(dataframe, column_name):
    """

    :param dataframe:
    :param column_name:
    :type column_name:str
    :return:
    """
    fig = plt.figure(figsize=(11.69, 8.27))
    p = plt.plot(dataframe.index, dataframe[column_name], 'b-', label=r"%s" % column_name)
    plt.hlines(0, min(dataframe.index), max(dataframe.index), 'r')
    plt.legend(loc='best')
    fig.autofmt_xdate(rotation=90)
    return p 

def available(show=True):
    if not show:
        return palettes.keys()
    else:
        f, ax = plt.subplots(5, 2, figsize=(5, 8))
        for i, name in enumerate(palettes.keys()):
            x, y = i // 2, i % 2
            cycle = palettes[name]
            for j, c in enumerate(cycle):
                ax[x, y].hlines(j, 0, 1, colors=c, linewidth=15)
            ax[x, y].set_ylim(-1, len(cycle))
            ax[x, y].set_title(name)
            despine(ax[x, y], True)
        plt.show() 

def Hlines(ys, x1, x2, **options):
    """Plots a set of horizontal lines.

    Args:
      ys: sequence of y values
      x1: sequence of x values
      x2: sequence of x values
      options: keyword args passed to plt.vlines
    """
    options = _UnderrideColor(options)
    options = _Underride(options, linewidth=1, alpha=0.5)
    plt.hlines(ys, x1, x2, **options) 

def example_plot_bias_ts():
    ax = pdf.plot(figsize=(8, 4), secondary_y='bias')
    plt.hlines(0, 1800, 2015, linestyles='-')
    ax.set_ylabel(r'$\mu$ (mm yr$^{-1}$ K$^{-1}$)')
    ax.set_title(r'$\mu$ candidates HEF')
    plt.ylabel(r'bias (mm yr$^{-1}$)')
    yl = plt.gca().get_ylim()
    plt.plot((res['t_star'], res['t_star']), (yl[0], 0),
             linestyle=':', color='grey')
    plt.ylim(yl)
    plt.tight_layout()
    plt.show() 

def test_limiter(self, default_calving):

        _, ds1, df_diag1 = default_calving

        model = FluxBasedModel(bu_tidewater_bed(),
                               mb_model=ScalarMassBalance(),
                               is_tidewater=True, calving_use_limiter=False,
                               flux_gate=0.06, do_kcalving=True,
                               calving_k=0.2)
        _, ds2 = model.run_until_and_store(3000)
        df_diag2 = model.get_diagnostics()
        assert_allclose(model.volume_m3 + model.calving_m3_since_y0,
                        model.flux_gate_m3_since_y0)
        assert_allclose(ds2.calving_m3[-1], model.calving_m3_since_y0)
        assert_allclose(ds2.volume_bsl_m3[-1], model.volume_bsl_km3 * 1e9)
        assert_allclose(ds2.volume_bwl_m3[-1], model.volume_bwl_km3 * 1e9)

        # Not exact same of course
        assert_allclose(ds1.volume_m3[-1], ds2.volume_m3[-1], rtol=0.06)
        assert_allclose(ds1.calving_m3[-1], ds2.calving_m3[-1], rtol=0.15)
        assert_allclose(ds1.volume_bsl_m3[-1], ds2.volume_bsl_m3[-1], rtol=0.3)
        assert_allclose(ds2.volume_bsl_m3, ds2.volume_bwl_m3)

        if do_plot:
            f, ax = plt.subplots(1, 1, figsize=(12, 5))
            df_diag1[['surface_h']].plot(ax=ax, color=['C3'])
            df_diag2[['surface_h', 'bed_h']].plot(ax=ax, color=['C1', 'k'])
            plt.hlines(0, 0, 60000, color='C0', linestyles=':')
            plt.ylim(-350, 800)
            plt.ylabel('Altitude [m]')
            plt.show() 

def show_iter_hist(fname, maxiter, nruns):
    """
    Helper routine to visualize the maximal iteration number across the simulation in a histogram

    Args:
        stats (dict): statistics object
        fname (str): filename
        maxiter: maximal iterations per run
        nruns: number of runs
    """

    # create plot and save
    fig, ax = plt.subplots(figsize=(15, 10))

    plt.hist(maxiter, bins=np.arange(min(maxiter), max(maxiter) + 2, 1), align='left', rwidth=0.9)

    # with correction allowed: axis instead of xticks
    # plt.axis([12, 51, 0, nruns+1])
    plt.xticks([13, 15, 20, 25, 30, 35, 40, 45, 50])

    ax.set_xlabel('iterations until convergence')

    plt.hlines(nruns, min(maxiter), max(maxiter), colors='red', linestyle='dashed')

    # with correction allowed: no logscale
    plt.yscale('log')

    plt.savefig(fname)

    assert os.path.isfile(fname), 'ERROR: plotting did not create PNG file' 

def test_water_level(self, default_calving):

        _, ds_1, _ = default_calving

        model = FluxBasedModel(bu_tidewater_bed(water_level=1000),
                               mb_model=ScalarMassBalance(),
                               is_tidewater=True, calving_use_limiter=True,
                               flux_gate=0.06, do_kcalving=True,
                               water_level=1000,
                               calving_k=0.2)
        _, ds_2 = model.run_until_and_store(3000)
        assert_allclose(model.volume_m3 + model.calving_m3_since_y0,
                        model.flux_gate_m3_since_y0)

        assert_allclose(ds_1.calving_m3, ds_2.calving_m3)

        if do_plot:
            df_diag = model.get_diagnostics()
            f, ax = plt.subplots(1, 1, figsize=(12, 5))
            df_diag[['surface_h', 'bed_h']].plot(ax=ax, color=['C3', 'k'])
            plt.hlines(1000, 0, 60000, color='C0', linestyles=':')
            plt.ylim(1000-350, 1000+800)
            plt.ylabel('Altitude [m]')
            plt.show()

        # Let the model decide the water level
        fls = bu_tidewater_bed(water_level=1000)
        thick = fls[0].thick
        thick[fls[0].bed_h > 1000] = cfg.PARAMS['free_board_lake_terminating']
        fls[0].thick = thick
        model = FluxBasedModel(fls, mb_model=ScalarMassBalance(),
                               is_tidewater=True, calving_use_limiter=True,
                               is_lake_terminating=True,
                               flux_gate=0.06, do_kcalving=True,
                               calving_k=0.2)
        assert_allclose(model.water_level, 1000, atol=1)

        with pytest.raises(InvalidParamsError):
            fls = bu_tidewater_bed(water_level=1000)
            FluxBasedModel(fls, mb_model=ScalarMassBalance(),
                           is_tidewater=True, calving_use_limiter=True,
                           is_lake_terminating=True) 

def plot_autocorrelation(data, plot=True, poly_order=1):
    """
    Given a dataset of a stationary vehicle on the ground,
    this compute the autocorrellation. The intersection with
    the horizontal line at 0.348 represents the correllation
    time constant. If the intersection does not occur, it
    indicates that the random walk process is not significant
    and is over powered by wide-band noise.
    """
    data.index = pd.TimedeltaIndex(data.index, unit='s')
    data_vals = []
    dt_vals = []
    dt = 1
    # downsample every dt seconds, taking mean
    data = data.resample(
        '{:d}L'.format(int(dt * 1000))).agg('mean')
    lag_max = int(dt * len(data.index) / 2)
    for i in range(1, int(lag_max / dt)):
        data_vals += [data.autocorr(i)]
        dt_vals += [i * dt]

    # polynomial fits
    p = np.polynomial.Polynomial.fit(dt_vals, data_vals, poly_order)

    # normalize by max of fit polynomial
    x = np.linspace(0, lag_max)
    y = p(x)
    y_max = np.max(np.abs(y))
    data_vals /= y_max
    y /= y_max
    p /= y_max

    if plot:
        plt.title('normalized autocorrelation')
        plt.plot(dt_vals[1:], data_vals[1:], '.', alpha=0.2, label='raw')
        plt.xlabel('lag, sec')
        plt.ylabel('corr/ max(corr)')
        plt.plot(x, y, linewidth=2, label='fit')
        plt.hlines(1 / np.e, 0, lag_max, linewidth=2)
        plt.gca().set_ylim(-1, 1)
        plt.gca().set_xlim(0, lag_max)
        plt.grid(True)

    correlation_time = _smallest_positive_real_root((p - 1 / np.e).roots())
    return correlation_time 

def visualize_outlierscore(value,label,contamination,path=None):
    """
    Visualize the predicted outlier score.

    Parameters
    ----------
    value: numpy array of shape (n_test, )
        The outlier score of the test data.
    label: numpy array of shape (n_test, )
        The label of test data produced by the algorithm.
    contamination : float in (0., 0.5), optional (default=0.1)
        The amount of contamination of the data set,
        i.e. the proportion of outliers in the data set. Used when fitting to
        define the threshold on the decision function.
    path: string
        The saving path for result figures.
    """

    sns.set(style="darkgrid")

    ts = np.arange(len(value))
    outlier_label=[]
    for i in range(len(ts)):
        if label[i]==1:
            outlier_label.append('inlier')
        else:
            outlier_label.append('outlier')
    X_outlier = pd.DataFrame({'ts':ts,'Outlier_score':value,'outlier_label':np.array(outlier_label)})
    pal = dict(inlier="#4CB391", outlier="gray")
    g = sns.FacetGrid(X_outlier, hue="outlier_label", palette=pal, height=5)
    g.map(plt.scatter, "ts", "Outlier_score", s=30, alpha=.7, linewidth=.5, edgecolor="white")

    ranking = np.sort(value)
    threshold = ranking[int((1 - contamination) * len(ranking))]
    plt.hlines(threshold, xmin=0, xmax=len(X_outlier)-1, colors="g", zorder=100, label='Threshold')
    threshold = ranking[int((contamination) * len(ranking))]
    plt.hlines(threshold, xmin=0, xmax=len(X_outlier)-1, colors="g", zorder=100, label='Threshold2')
    if path:
        plt.savefig(path+'/visualize_outlierscore.png')
    plt.show()



# def visualize_outlierresult(X,label,path=None):
#     """
#     Visualize the predicted outlier result.
#
#     Parameters
#     ----------
#     X: numpy array of shape (n_test, n_features)
#         The test data.
#     label: numpy array of shape (n_test, )
#         The label of test data produced by the algorithm.
#
#     """
#     X['outlier']=pd.Series(label)
#     pal = dict(inlier="#4CB391", outlier="gray")
#     g = sns.pairplot(X, hue="outlier", palette=pal)
#     if path:
#         plt.savefig(path+'/visualize_outlierresult.png')
#     plt.show() 

def _draw_figs(cls, topo_tensors, alloc_plan, cmap, figsize, fontsize, lw, split_on_large_graph, num_tensors_per_split):
    xmins = [alloc_plan.plan[tensor_name].offset_start for tensor_name in topo_tensors]
    xmaxs = [alloc_plan.plan[tensor_name].offset_end for tensor_name in topo_tensors]
    colors = [cmap(random()) for _ in alloc_plan.plan]
    labels = topo_tensors[:]
    sizes = [alloc_plan.plan[tensor_name].size for tensor_name in topo_tensors]
    if split_on_large_graph:
      xmin_chunks = [xmins[i:i+num_tensors_per_split] for i in range(0, len(xmins), num_tensors_per_split)]
      xmax_chunks = [xmaxs[i:i+num_tensors_per_split] for i in range(0, len(xmaxs), num_tensors_per_split)]
      color_chunks = [colors[i:i+num_tensors_per_split] for i in range(0, len(colors), num_tensors_per_split)]
      label_chunks = [labels[i:i+num_tensors_per_split] for i in range(0, len(labels), num_tensors_per_split)]
      size_chunks = [sizes[i:i+num_tensors_per_split] for i in range(0, len(sizes), num_tensors_per_split)]
    else:
      xmin_chunks = [xmins]
      xmax_chunks = [xmaxs]
      color_chunks = [colors]
      label_chunks = [labels]
      size_chunks = [sizes]
    figs = []
    for i, (xmins, xmaxs, colors, labels, sizes) in enumerate(zip(xmin_chunks, xmax_chunks, color_chunks, label_chunks, size_chunks)):
      fig, _ = plt.subplots(1, 1)
      ys = [len(xmins)-i for i in range(len(xmins))]
      for y, xmin, xmax, color, size in zip(ys, xmins, xmaxs, colors, sizes):
        plt.hlines(y, xmin, xmax, lw=lw, colors=color)
        plt.text(xmax+lw*10, y-0.01*lw, '{} bytes'.format(size), fontdict={'fontsize': fontsize})
      plt.xlabel('Offset (bytes)', fontdict={'fontsize': fontsize})
      plt.yticks(ys, labels, fontsize=fontsize)
      plt.xticks(fontsize=fontsize)
      plt.ylabel(
        'Tensor Names (Topological Ordered, Top to Bottom)',
        fontdict={'fontsize':fontsize}
      )
      if i:
        title = 'Optimal Tensor Allocation: {} bytes in total (Cont.)'.format(alloc_plan.total_size)
      else:
        title = 'Optimal Tensor Allocation: {} bytes in total'.format(alloc_plan.total_size)
      plt.title(
        title,
        fontdict={'fontsize': fontsize}
      )
      if figsize is None:
        figsize = (num_tensors_per_split, num_tensors_per_split / 2)
      fig.set_size_inches(*figsize)
      fig.tight_layout()
      figs.append(fig)
    return figs 

def plot_best_state(train_state):
    X_train, y_train, X_test, y_test, best_y, best_ystd = train_state.packed_data()

    y_dim = best_y.shape[1]

    # Regression plot
    if X_test.shape[1] == 1:
        idx = np.argsort(X_test[:, 0])
        X = X_test[idx, 0]
        plt.figure()
        for i in range(y_dim):
            if y_train is not None:
                plt.plot(X_train[:, 0], y_train[:, i], "ok", label="Train")
            if y_test is not None:
                plt.plot(X, y_test[idx, i], "-k", label="True")
            plt.plot(X, best_y[idx, i], "--r", label="Prediction")
            if best_ystd is not None:
                plt.plot(
                    X, best_y[idx, i] + 2 * best_ystd[idx, i], "-b", label="95% CI"
                )
                plt.plot(X, best_y[idx, i] - 2 * best_ystd[idx, i], "-b")
        plt.xlabel("x")
        plt.ylabel("y")
        plt.legend()
    elif X_test.shape[1] == 2:
        for i in range(y_dim):
            plt.figure()
            ax = plt.axes(projection=Axes3D.name)
            ax.plot3D(X_test[:, 0], X_test[:, 1], best_y[:, i], ".")
            ax.set_xlabel("$x_1$")
            ax.set_ylabel("$x_2$")
            ax.set_zlabel("$y_{}$".format(i + 1))

    # Residual plot
    if y_test is not None:
        plt.figure()
        residual = y_test[:, 0] - best_y[:, 0]
        plt.plot(best_y[:, 0], residual, "o", zorder=1)
        plt.hlines(0, plt.xlim()[0], plt.xlim()[1], linestyles="dashed", zorder=2)
        plt.xlabel("Predicted")
        plt.ylabel("Residual = Observed - Predicted")
        plt.tight_layout()

    if best_ystd is not None:
        plt.figure()
        for i in range(y_dim):
            plt.plot(X_test[:, 0], best_ystd[:, i], "-b")
            plt.plot(
                X_train[:, 0],
                np.interp(X_train[:, 0], X_test[:, 0], best_ystd[:, i]),
                "ok",
            )
        plt.xlabel("x")
        plt.ylabel("std(y)") 

def _plot(self, series, forecasts, n_test, file_name, sensor_name, time, n_seq):
        """
        Same as function 'plot_forecasts', replace the datetime in x-axis with days.
        """
        plot_one_line = 1
        label_fontsize = 35
        axis_fontsize = 30
        linewidth = 5

        # plot the entire dataset in blue
        fig = pyplot.figure()
        ax1 = fig.add_subplot(1, 1, 1)
        # make x label in a specific format
        # ax1.xaxis_date()
        # ax1.xaxis.set_major_formatter(DateFormatter('%m-%d'))
        forecasts = np.array(forecasts)
        pyplot.plot(series.index, series.values, label='Actual data', linewidth=linewidth)
        ####################### plot the forecast value #########################
        X = []
        for i in range(1, forecasts.shape[1] + 1):
            off_s = len(series) - n_test + i - n_seq
            off_e = off_s + n_test - 1
            X.append(range(off_s, off_e + 1))
        X = np.array(X)
        Y = np.array(forecasts)
        for i in range(0, Y.shape[1]):
            index = X[i]
            pyplot.plot(np.arange(index[0], index[-1] + 1), Y[:, i], label='Prediction: t+' + str(i + 1),
                        linewidth=linewidth)
            if plot_one_line == 1:
                break
        pyplot.hlines(self.threshold[self.sensor_name][0], series.index[0], series.index[-1], colors='r',
                      linewidth=linewidth)
        pyplot.hlines(self.threshold[self.sensor_name][1], series.index[0], series.index[-1], colors='g', label='normal',
                      linewidth=linewidth)
        pyplot.hlines(self.threshold[self.sensor_name][2], series.index[0], series.index[-1], colors='r',
                      linewidth=linewidth)

        pyplot.title(self.sensor_name_acronym[self.sensor_name], fontsize=label_fontsize)
        pyplot.legend(fontsize=label_fontsize, loc='upper left')
        pyplot.xlabel('Days', fontsize=label_fontsize)
        pyplot.ylabel(self.units[sensor_name], fontsize=label_fontsize)
        pyplot.xticks(fontsize=axis_fontsize)
        pyplot.yticks(fontsize=axis_fontsize)
        # replace date to sequential days

        # show the plot
        fig.show()

        if self.save_info:
            fig.set_size_inches(18.5, 10.5)
            fig.savefig(os.path.join(self.file_path, file_name + '.png'), bbox_inches='tight', dpi=150)

        pyplot.close(fig)