def plot_2d_clusters(X, labels, centers):
    """
    Given an observation array, a label vector, and the location of the centers
    plot the clusters
    """

    clabels = set(labels)
    K = len(clabels)

    if len(centers) != K:
        raise ValueError("Expecting the number of unique labels and centres to"
                         " be the same!")

    # Plot the true clusters
    figure(figsize=(10, 10))
    ax = gca()

    vor = Voronoi(centers)

    voronoi_plot_2d(vor, ax)

    colors = cm.hsv(np.arange(K)/float(K))
    for k, col in enumerate(colors):
        my_members = labels == k
        scatter(X[my_members, 0], X[my_members, 1], c=col, marker='o', s=20)

    for k, col in enumerate(colors):
        cluster_center = centers[k]
        scatter(cluster_center[0], cluster_center[1], c=col, marker='o', s=200)

    axis('tight')
    axis('equal')
    title('Clusters')

def pathsByLength(geo, targ_dict, background=True):
  """
  Overlay paths all on one tree colored by path length.
  """
  lengths = [targ_dict[s]['pathLength'] for s in targ_dict.keys()]
  plt.figure(figsize=(4,6))
  
  if background:
    branchpts = []
    # Plot by branches
    for b in geo.branches: 
      branchpts.append([[n.x, n.y] for n in b.nodes])
    for b in branchpts: # Plot the background skeleton first
        for s in range(len(b)-1):
          plt.plot([b[s][0], b[s+1][0]],
                   [b[s][1], b[s+1][1]], color='black', alpha=0.3)
  
  # Plot each thing
  pDF = PathDistanceFinder(geo, geo.soma)
  for t in targ_dict.keys():
    targ = [seg for seg in geo.segments if 
              seg.filamentIndex==targ_dict[t]['closestFil']][0]
    path = pDF.pathTo(targ)    # This makes sure the colors range the whole spectrm
    pcolor = (targ_dict[t]['pathLength']-min(lengths))/(max(lengths)-min(lengths))
    for p in range(len(path)-1):
      c0, c1 = path[p].coordAt(0), path[p+1].coordAt(0)
      plt.plot([c0[0], c1[0]], [c0[1], c1[1]], linewidth=3,
               color=cm.hsv(pcolor), alpha=0.7)
  plt.show()

def histogram(names, values, title="", xlabel="", ylabel=""):

    if not names or not values:
        return _empty_chart(title, xlabel, ylabel)

    figure = mpl_Figure()
    canvas = mpl_FigureCanvas(figure)
    figure.set_canvas(canvas)

    ax = figure.add_subplot(111, projection=BasicChart.name)

    colors = [cm.hsv(float(i)/len(values)) for i in xrange(len(values))]
    n, bins, patches = ax.hist(
        values, 10, normed=0, histtype="bar", label=names, color=colors)

    for label in ax.xaxis.get_ticklabels():
        label.set_rotation(-35)
        label.set_horizontalalignment('left')

    ax.plegend = ax.legend(loc="upper right", fancybox=True, shadow=True)

    ax.xaxis.set_major_formatter(mpl_FuncFormatter(
        lambda time, pos: utils.time_to_string(time)[:-7]))
    ax.set_xlim(xmin=0)
    ax.set_title(title)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)

    return ChartWidget(figure, xlock=True)

def plot_spectra_grid(file_set,protein,ligands,ligand):
    grid = len(protein) + len(ligand)
    
    # pick the correct file
    proteins = file_set.keys()
    index = ligands.index(ligand)
    file = file_set[protein][index]
    
    # pick a title
    title = "%s - %s" %(protein, ligand)
    
    # make a dataframe
    df = xml2df(file)
    
    # plot the spectra
    fig = plt.figure();
    ax = df['fluorescence'].iloc[:,12].plot(ylim=(0,100000),legend=False, linewidth=4,color='m');
    ax.axvline(x=480,color='0.7',linestyle='--');
    for i in range(11):
        #s = df['fluorescence'].iloc[:,i].plot(ylim=(0,100000),linewidth=3,c=cm.hsv(i*15), ax = ax, title=title);
        df['fluorescence'].iloc[:,i].plot(ylim=(0,100000),linewidth=3,c=cm.hsv(i*15), ax = ax);
        df['fluorescence'].iloc[:,11+i].plot(ylim=(0,100000),legend=False, linewidth=4,c=cm.gray(i*15+50), ax = ax, fontsize =20);
    sns.despine()
    plt.xlim(320,600)
    plt.yticks([])
    plt.xlabel('wavelength (nm)', fontsize=20)
    plt.tight_layout();
    plt.savefig('%s.eps'%title, type='eps', dpi=1000)

def animate():
    global graph, pfad, kreis, ctr, cnt

    # neuen Punkt generieren
    pts, R, attached = dla()

    # Pfad und Kreis updaten
    x,y = zip(*pts)
    pfad.set_data(x, y)
    kreis.set_data([ctr + R*cos(phi) for phi in linspace(0,2*pi,100)],
                   [ctr + R*sin(phi) for phi in linspace(0,2*pi,100)])

    # und letzten Punkt permanent ausgeben, falls angelagert
    if attached:
        graph.scatter(x[-1], y[-1], s=4, c=cm.hsv(cnt), linewidth=0)
        cnt += 0.1

    # Groesse anpasse
    graph.axis((0,M-1,0,M-1))

    # periodisch die Funktion wieder aufrufen, wenn noch ok
    if R >= 0.4*M:
        print "Radius %d ist zu gross geworden" % R
    else:
        figure.canvas.manager.window.after(200, animate)

    figure.canvas.draw()

def command(paths):
    """ Show a lowest elevation image. """
    # order
    index = dict(
        NL60=0,
        NL61=1,
        nhb=2,
        ess=3,
        emd=4,
        JAB=5,
    )

    # load
    d_rang, d_elev, d_anth = {}, {}, {}
    for p in paths:
        with h5py.File(p, 'r') as h5:
            for r in h5:
                if r in d_rang or r == 'ase':
                    continue
                d_rang[r] = h5[r]['range'][:]
                d_elev[r] = h5[r]['elevation'][:]
                d_anth[r] = h5[r].attrs['antenna_height']
    radars = d_anth.keys()
    elev = np.ma.empty((len(radars),) + d_rang[radars[0]].shape)
    rang = np.ma.empty((len(radars),) + d_rang[radars[0]].shape)
    anth = np.empty((len(radars), 1, 1))
    for r in radars:
        elev[index[r]] = np.ma.masked_equal(d_elev[r], config.NODATAVALUE)
        rang[index[r]] = np.ma.masked_equal(d_rang[r], config.NODATAVALUE)
        anth[index[r]] = float(d_anth[r]) / 1000

    # calculate
    theta = calc.calculate_theta(
        rang=rang, elev=np.radians(elev), anth=anth,
    )
    alt = calc.calculate_height(
        theta=theta, elev=np.radians(elev), anth=anth,
    )
    which = np.ma.where(
        alt == alt.min(0),
        np.indices(alt.shape)[0],
        -1,
    ).max(0)
    what = alt.min(0)

    # colors
    hue = cm.hsv(colors.Normalize(vmax=len(radars))(which), bytes=True)
    sat = 1 - colors.Normalize(vmax=5, clip=True)(what)[..., np.newaxis]

    hue[..., :3] *= sat
    hue[sat.mask[..., 0]] = 255
    Image.fromarray(hue).save('elevation_image.png')
    return 0

def create_pred_movie_no_conf(conf, predList, moviename, outmovie, outtype):
    predLocs, predscores, predmaxscores = predList
    #     assert false, 'stop here'
    tdir = tempfile.mkdtemp()

    cap = cv2.VideoCapture(moviename)
    nframes = int(cap.get(cvc.FRAME_COUNT))

    cmap = cm.get_cmap('jet')
    rgba = cmap(np.linspace(0, 1, conf.n_classes))

    fig = mpl.figure.Figure(figsize=(9, 4))
    canvas = FigureCanvasAgg(fig)
    for curl in range(nframes):
        framein = myutils.readframe(cap, curl)
        framein = crop_images(framein, conf)

        fig.clf()
        ax1 = fig.add_subplot(1, 2, 1)
        ax1.imshow(framein[:, :, 0], cmap=cm.gray)
        ax1.scatter(predLocs[curl, :, 0, 0], predLocs[curl, :, 0, 1],  # hold=True,
                    c=cm.hsv(np.linspace(0, 1 - old_div(1., conf.n_classes), conf.n_classes)),
                    s=20, linewidths=0, edgecolors='face')
        ax1.axis('off')
        ax2 = fig.add_subplot(1, 2, 2)
        if outtype == 1:
            curpreds = predscores[curl, :, :, :, 0]
        elif outtype == 2:
            curpreds = predscores[curl, :, :, :, 0] * 2 - 1

        rgbim = create_pred_image(curpreds, conf.n_classes)
        ax2.imshow(rgbim)
        ax2.axis('off')

        fname = "test_{:06d}.png".format(curl)

        # to printout without X.
        # From: http://www.dalkescientific.com/writings/diary/archive/2005/04/23/matplotlib_without_gui.html
        # The size * the dpi gives the final image size
        #   a4"x4" image * 80 dpi ==> 320x320 pixel image
        canvas.print_figure(os.path.join(tdir, fname), dpi=80)

        # below is the easy way.
    #         plt.savefig(os.path.join(tdir,fname))

    tfilestr = os.path.join(tdir, 'test_*.png')
    mencoder_cmd = "mencoder mf://" + tfilestr + " -frames " + "{:d}".format(
        nframes) + " -mf type=png:fps=15 -o " + outmovie + " -ovc lavc -lavcopts vcodec=mpeg4:vbitrate=2000000"
    os.system(mencoder_cmd)
    cap.release()

    def _on_update(self, histogram, rois):
        if histogram not in self.map:
            return
        source, s_out = self.map[histogram]
        N = len(rois)
        data = source()

        gray = self._make_gray(data)
        s_out.data = regular_array2rgbx(gray)
        for i in range(N):
            color = (255 * plt_cm.hsv([float(i) / max(N, 10)])).astype('uint8')
            color = regular_array2rgbx(color)
            r = rois[i]
            mask = (data < r.right) & (data >= r.left)
            s_out.data[mask] = color
        s_out.update_plot()

def plot_mse(mse_list,label_list):
    for i,mse in enumerate(mse_list):
        c = cm.hsv(np.linspace(0,1,len(mse_list)+1))
        data_x = range(mse.shape[0])
        data_y = mse;
        plt.plot(data_x,data_y,color=c[i,:])
    plt.title('Sparse autoencoder training curve MSE')
    plt.xlabel('Epoch')
    plt.ylabel('MSE')
    plt.legend(label_list,'upper right')
    plt.grid(True)
    plt.xlim(0,2000)
    plt.ylim(0,100)
    f = plt.gcf();
    f.set_size_inches(19.2,10.8)
    plt.savefig('../result_images/mnist_dictionary_training.png',dpi=100)
    plt.close()

def plot_seimicity_rate(earthquakes, time = 'hour', figsize = (8,8)):
    '''
    Function get from Thomas Lecocq
    http://www.geophysique.be/2013/09/25/seismicity-rate-using-obspy-and-pandas/
    '''
    
    m_range = (-1,11)
    time = time.lower()
    
    if time == 'second':
        time_format = '%y-%m-%d, %H:%M:%S %p'
    elif time == 'minute':
        time_format = '%y-%m-%d, %H:%M %p'
    elif time == 'hour':
        time_format = '%y-%m-%d, %H %p'
    elif time == 'day':
        time_format = '%y-%m-%d'
    elif time == 'month':
        time_format = '%y-%m'
    elif time == 'year':
        time_format = '%y'
    else:
        time_format = '%y-%m-%d, %H %p'
    
    df = eq2df(earthquakes)

    #Arrange quakes by their magnitude and color appropiately
    
    bins = np.arange(m_range[0], m_range[1])
    labels = np.array(["[%i:%i)"%(i,i+1) for i in bins])
    colors = [cm.hsv(float(i+1)/(len(bins)-1)) for i in bins]
    df['Magnitude_Range'] = pd.cut(df['mag'], bins=bins, labels=False)
    
    df['Magnitude_Range'] = labels[df['Magnitude_Range'].values]
    
    df['Year_Month_day'] = [di.strftime(time_format) for di in df.index]
    
    rate = pd.crosstab(df.Year_Month_day, df.Magnitude_Range)
    
    rate.plot(kind='bar',stacked=True,color=colors,figsize=figsize)
    plt.legend(loc='best')
    plt.ylabel('Number of earthquakes')
    plt.xlabel('Date and Time')
    plt.tight_layout()
    plt.grid()
    plt.show()

def gen_colors(clusters):
    """
    Takes 'clusters' and creates a list of colors so a cluster has a color.
    
    :param clusters: A flat list of integers where an integer represents which
        cluster the information belongs to.
    :type clusters: list
    
    :returns: colorm : list
        A list of colors obtained from matplotlib colormap cm.hsv. The
        length of 'colorm' is the same as the number of distinct
        clusters.
    """
    import matplotlib.cm as cm
    
    n = len(set(clusters))
    colorm = [cm.hsv(i * 1.0 /n, 1) for i in xrange(n)]
    return colorm

    def DisplaySavedMapFrame(self, worldFrame, resourcesGRMaxE, frameNo, colours, creatSpec):
        TileBlitSize = (int(self.tw*(self.RESIZE[0]/float(self.SIZE[0]))), int(self.tw*(self.RESIZE[1]/float(self.SIZE[1]))))
        sizeRat0 = self.RESIZE[0]/float(self.SIZE[0])
        sizeRat1 = self.RESIZE[1]/float(self.SIZE[1])
        done = False
        step = 0
        while not done:
            for event in pygame.event.get():
                if event.type == pygame.QUIT:
                    done = True
                    pygame.quit()
                    return
            if step == 0:
                t0 = time.time()
                self.screen.fill(self.WHITE)
                resources = worldFrame[2]
                for a,b in np.ndindex((resources.shape[0], resources.shape[1])):
                    x, y = self.gridWidth-a-1, self.gridHeight-b-1
                    image = pygame.transform.scale(self.gameMap.get_tile_image(x,y,0), TileBlitSize) #pre-scaled at load time
                    self.screen.blit(image, self.TransformMap(np.array([x, y])))
                    if resourcesGRMaxE[0][x][y] > 0:
                        if len(worldFrame) == 4:
                            resourcesGRMaxE[1][x][y] *= worldFrame[3]
                        tempTransPos = self.TransformResourcePos(np.array([x,y]))
                        polygonPos = [(tempTransPos[0], tempTransPos[1]+int((self.th/16.)*sizeRat1)), (tempTransPos[0]+int((self.tw*7./16.)*sizeRat0), tempTransPos[1]+int((self.th/2.)*sizeRat1)), (tempTransPos[0], tempTransPos[1]+int((self.th*15./16.)*sizeRat1)), (tempTransPos[0]-int((self.tw*7./16.)*sizeRat0), tempTransPos[1]+int((self.tw/4.)*sizeRat1))]
                        polygonColour = self.RED[0]*worldFrame[2][x][y]/float(resourcesGRMaxE[1][x][y])
                        pygame.draw.aalines(self.screen, [polygonColour, 0, 0], True, polygonPos, 1)
                        #pygame.draw.polygon(self.screen, PolygonColourArray[x][y][step], PolygonPosArray[x][y], 1)
                t1 = time.time()

                for i in xrange(len(worldFrame[0])):
                    for j in xrange(len(creatSpec)):
                        if worldFrame[0][i][0] == creatSpec[j][0]:
                            creaturePos = self.TransformPos(np.array([worldFrame[0][i][1], worldFrame[0][i][2]]))
                            creatureColour = cm.hsv(colours[j], bytes=True)[0:3]
                            pygame.draw.circle(self.screen, self.WHITE, creaturePos, 5)
                            pygame.draw.circle(self.screen, self.BLACK, creaturePos, 4)
                            pygame.draw.circle(self.screen, creatureColour, creaturePos, 2)
                t2 = time.time()
                
                print('Frame time = %s, for1 = %s, for2 = %s' % (time.time()-t0, t1-t0, t2-t1))
                pygame.display.flip()
            self.clock.tick(15)
            step += 1

def plot_G(ax, Gs, weights, intersect):
    
    G_x=np.arange(0,5,.001)
    l = len(Gs)
    G_ys = []
    for i in xrange(l):
        c = cm.hsv(float(i)/l,1)
        mu = Gs[i][0]
        var = Gs[i][1]
        
        G_y = mlab.normpdf(G_x, mu, var**.5)*weights[i]
        G_ys.append(G_y)
        ax.plot(G_x,G_y,color=c)
        ax.plot(mu,-.001,"^",ms=10,alpha=.7,color=c)
    
    #ax.plot(G_x,np.power(G_ys[1]-G_ys[0],1),color='k')
    if intersect[0] !=None:
        ax.plot(intersect[0],intersect[1],"|",ms=10,alpha=.7,color='k')
    ax.plot([0,5],[0,0],color='k')

def plot_seimicity_rate(earthquakes, time="hour", figsize=(12, 8)):
    """
    Function get from Thomas Lecocq
    http://www.geophysique.be/2013/09/25/seismicity-rate-using-obspy-and-pandas/
    """

    m_range = (-1, 11)
    time = time.lower()

    if time == "second":
        time_format = "%y-%m-%d, %H:%M:%S %p"
    elif time == "minute":
        time_format = "%y-%m-%d, %H:%M %p"
    elif time == "hour":
        time_format = "%y-%m-%d, %H %p"
    elif time == "day":
        time_format = "%y-%m-%d"
    elif time == "month":
        time_format = "%y-%m"
    elif time == "year":
        time_format = "%y"
    else:
        time_format = "%y-%m-%d, %H %p"

    df = eq2df(earthquakes)

    bins = np.arange(m_range[0], m_range[1])
    labels = np.array(["[%i:%i)" % (i, i + 1) for i in bins])
    colors = [cm.hsv(float(i + 1) / (len(bins) - 1)) for i in bins]
    df["Magnitude_Range"] = pd.cut(df["mag"], bins=bins, labels=False)

    df["Magnitude_Range"] = labels[df["Magnitude_Range"].values]

    df["Year_Month_day"] = [di.strftime(time_format) for di in df.index]

    rate = pd.crosstab(df.Year_Month_day, df.Magnitude_Range)

    rate.plot(kind="bar", stacked=True, color=colors, figsize=figsize)
    plt.legend(bbox_to_anchor=(1.20, 1.05))
    plt.ylabel("Number of earthquakes")
    plt.xlabel("Date and Time")
    plt.show()

    def simple_plot(self, fn_out):

        cps = self.mus
        plt.rc('grid',color='0.75',linestyle='l',linewidth='0.1')
        fig, ax_arr = plt.subplots(1,3)
        fig.set_figwidth(9)
        fig.set_figheight(4)
        axescolor  = '#f6f6f6'
        print ax_arr 
        print self.bics
        print self.params

        ax_arr[1].plot(self.params, self.bics)

        n, bins, patches = ax_arr[0].hist(cps,alpha=.9,ec='none',normed=1,color='#8DABFC',bins=len(cps)/10)
        #self.addGMM(gX.gmm, axarr[1,1], cps, gX.labels, overlaps)
        
        G_x=np.arange(0,max(cps)+1,.1)
        l = self.gmm.means.shape[0]
        
        for i in xrange(l):
            c = cm.hsv(float(i)/l,1)
            mu = self.gmm.means[i,0]
            #var = self.gmm.covars[i][0][0]
            var = self.gmm.covars[i]

            G_y = mlab.normpdf(G_x, mu, var**.5)*self.gmm.weights[i]
            ax_arr[0].plot(G_x,G_y,color=c)
            ax_arr[0].plot(mu,-.001,"^",ms=10,alpha=.7,color=c)
        
        if np.amax(cps)<2:
            ax_arr[0].set_xlim(0,2)
        ylims = ax_arr[0].get_ylim()
        if ylims[1] > 10:
            ax_arr[0].set_ylim(0,10)

        fig.sca(ax_arr[2]) 
        dendro = hclust.dendrogram(self.Z, orientation='right')
        ylims = ax_arr[2].get_ylim()
        ax_arr[2].set_ylim(ylims[0]-1, ylims[1]+1)

        fig.savefig(fn_out)

    def draw_data(self, data):
        totalsize = sum(zip(*data)[1])     # get sum of subentry sizes
        unit = 1.5 * np.pi / totalsize

        angle = 0.5 * np.pi
        if self.log:
            maxy = max(map(np.log2, zip(*data)[2]))
        else:
            maxy = max(zip(*data)[2])
        for d in data:
            relangle = unit * d[1]

            if len(d[3]) > 0 or self.mode == 'user':
                # scale colours since the legend occupies a quarter
                scaledangle = (angle - 0.5*np.pi) * (2 / 1.5)
                colour = cm.hsv(scaledangle/(2*np.pi))
            else:
#                colour = cm.Greys(scaledangle/(2*np.pi))
                colour = "#999999"

            if self.log:
                # take logarithm to accomodate for big differences
                y = np.log2(d[2])
            else:
                y = d[2]

            bar = ax.bar(angle, y, width=relangle, bottom=maxy*0.2,
                         color=colour, label=d[0],
                         picker=True)
            angle += relangle

            if self.mode == 'dir':
                desc = '{0}\n{1}G'.format(d[0], d[1])
            elif self.mode == 'user':
                desc = '{0}\n{1}%'.format(d[0], d[1])
            self.draw_desc(bar[0], d, desc)

        self.draw_legend(maxy)

        fig.canvas.draw()

def plot_2d_GMMs(X, labels, means, covs, percentcontour=0.66, npoints=30):
    """
    Given an observation array, a label vector (integer values), and GMM mean
    and covariance parameters, plot the clusters and parameters.
    """

    clabels = set(labels)
    K = len(clabels)

    if len(means) != len(covs) != K:
        raise ValueError("Expecting the number of unique labels, means and"
                         "covariances to be the same!")

    phi = np.linspace(-np.pi, np.pi, npoints)

    circle = np.array([np.sin(phi), np.cos(phi)]).T

    figure(figsize=(10, 10))
    gca()

    colors = cm.hsv(np.arange(K)/float(K))
    for k, col in zip(clabels, colors):

        # points
        my_members = labels == k
        scatter(X[my_members, 0], X[my_members, 1], c=col, marker='o', s=20)

        # means
        cluster_center = means[k, :]
        scatter(cluster_center[0], cluster_center[1], c=col, marker='o', s=200)

        # covariance
        L = la.cholesky(np.array(covs[k]) * chi2.ppf(percentcontour, [3])
                        + 1e-5 * np.eye(covs[k].shape[0]))
        covpoints = circle.dot(L) + means[k, :]
        plot(covpoints[:, 0], covpoints[:, 1], color=col, linewidth=3)

    axis('tight')
    axis('equal')
    title('Clusters')

def plot_distances(x, y, labels):
    print "----"
    print labels
    print x
    print y
    print "----"

    fig, ax = plt.subplots()
    colors = np.random.rand(21)  # 21 datasets
    dataset_dict = dict()
    markers = dict()
    markers['RAND'] = '>'
    markers['SMAC'] = 'o'
    markers['TPE'] = (5, 1)

    # shape by strategy 
    # color by dataset
    color_index = 0
    for p in range(0, len(x)):
        tmp = labels[p].split('.')
        dataset = tmp[0]

        dataset_index = datasets.index(dataset)
        color = cm.hsv(dataset_index / 25., 1)

        strategy = tmp[1]
        plt.plot(x[p], y[p], marker=markers[strategy], c=color, alpha=.75, label=labels[p])

    # plt.tight_layout()
    plt.title('error variance vs configuration dissimilarity')
    plt.xlabel('MCPS dissimilarity')
    plt.ylabel('Error variance')
    # TODO fix legend
    # handles, labels = ax.get_legend_handles_labels()
    # lgd = ax.legend(handles, labels, loc='center right', bbox_to_anchor=(0.5,-0.1), numpoints=1)
    plt.savefig('../distances%s/_all.png' % suffix, dpi=200, bbox_inches='tight')
    plt.close("all")

import numpy as np
import matplotlib.pylab as plt
import matplotlib.cm as cm

loci = 10000
runs = 200
data = np.zeros((runs,loci))

i = 0
for line in open('out.txt', 'rb'):
    hist = line.split()
    data[i,int(hist[0])-1] = int(hist[2])
    if (int(hist[0]) == loci):
        i+=1

for i in range(15):
    print sum(data[i,:])
    plt.plot(range(loci),data[i,:], 'o', color=cm.hsv(i/15.,1), label='Generation %i'%(5*i))
plt.xlabel("Block Width")
plt.ylabel("Count")
plt.legend()
plt.show()

def differential_displacement(base_system, disl_system, burgers_vector, plot_range, 
       neighbor_list = None, neighbor_list_cutoff = None, component = 'standard',
       axes = None, plot_scale = 1, save_file = None, show = True):
    """
    Function for generating a differential displacement plot for a 
    dislocation containing system.
    
    Arguments:
    base_system -- atomman.System defect-free reference system corresponding 
                   to disl_system.
    disl_system -- atomman.System system containing the defect.
    burgers_vector -- 3x1 numpy array for the dislocation's Burgers vector.
    plot_range -- 3x3 numpy array specifying the Cartesian space to include 
                  atoms in the plot.
    
    Optional Keyword Arguments:
    neighbor_list -- pre-computed neighbor list for base_system.
    neighbor_list_cutoff -- cutoff for computing a neighbor list for 
                            base_system.
    component -- indicates the style of the calculation to use.
    axes -- 3x3 numpy array indicating the crystallographic 
                             axes corresponding to the box's Cartesian axes. 
                             If given, only used for transforming the 
                             burgers_vector.
    plot_scale -- scalar for multiplying the magnitude of the differential 
                  displacement arrows.
    save_file -- if given then the plot will be saved to a file with this name.
    show -- Boolean flag for showing the figure. Default is True.  
    """
    #Burgers vector setup
    if axes is not None:
        T = am.tools.axes_check(axes)
        burgers_vector = T.dot(burgers_vector)
    burgers_vector_magnitude = np.linalg.norm(burgers_vector)
    burgers_vector_uvect = burgers_vector / burgers_vector_magnitude    
    
    #neighbor list setup
    if neighbor_list is not None:
        assert neighbor_list_cutoff is None, 'neighbor_list and neighbor_list_cutoff cannot both be given'
    elif neighbor_list_cutoff is not None:
        neighbor_list = am.NeighborList(base_system, neighbor_list_cutoff)
    elif 'neighbors' in base_system.prop:
        neighbor_list = base_system.prop['neighbors']
    elif 'nlist' in base_system.prop:
        neighbor_list = base_system.prop['nlist']
    
    if isinstance(neighbor_list, am.NeighborList):
        objectnlist = True
    else:
        objectnlist = False
    
    #Identify atoms in plot range
    base_pos = base_system.atoms_prop(key='pos')
    plot_range_indices = np.where((base_pos[:, 0] > plot_range[0,0]) & (base_pos[:, 0] < plot_range[0,1]) & 
                                  (base_pos[:, 1] > plot_range[1,0]) & (base_pos[:, 1] < plot_range[1,1]) &
                                  (base_pos[:, 2] > plot_range[2,0]) & (base_pos[:, 2] < plot_range[2,1]))[0]
    
    #initial plot setup and parameters
    fig, ax1, = plt.subplots(1, 1, squeeze=True, figsize=(7,7), dpi=72)
    ax1.axis([plot_range[0,0], plot_range[0,1], plot_range[1,0], plot_range[1,1]])
    atom_circle_radius = burgers_vector_magnitude / 10
    arrow_width_scale = 1. / 200.

    #Loop over all atoms i in plot range
    for i in plot_range_indices:

        #Plot a circle for atom i
        color = cm.hsv((base_pos[i, 2] - plot_range[2,0]) / (plot_range[2,1] - plot_range[2,0]))
        ax1.add_patch(mpatches.Circle(base_pos[i, :2], atom_circle_radius, fc=color, ec='k'))
    
        #make list of all neighbors for atom i
        if objectnlist:
            neighbor_indices = neighbor_list[i]
        else:
            neighbor_indices = neighbor_list[i, 1 : neighbor_list[i, 0] + 1]

        #Compute distance vectors between atom i and its neighbors for both systems        
        base_dvectors = base_system.dvect(int(i), neighbor_indices)
        disl_dvectors = disl_system.dvect(int(i), neighbor_indices)
    
        #Compute differential displacement vectors
        dd_vectors = disl_dvectors - base_dvectors
        
        #Compute centerpoint positions for the vectors
        arrow_centers = base_pos[i] + base_dvectors / 2
        
        if component == 'standard':
            #compute unit distance vectors
            base_uvectors = base_dvectors / np.linalg.norm(base_dvectors, axis=1)[:,np.newaxis]

            #compute component of the dd_vector parallel to the burgers vector
            dd_components = dd_vectors.dot(burgers_vector_uvect)        
            dd_components[dd_components > burgers_vector_magnitude / 2] -= burgers_vector_magnitude
            dd_components[dd_components < -burgers_vector_magnitude / 2] += burgers_vector_magnitude
            
            #scale arrow lengths and vectors
            arrow_lengths = base_uvectors * dd_components[:,np.newaxis] * plot_scale
            arrow_widths = arrow_width_scale * dd_components * plot_scale
        
            #plot the arrows
#.........这里部分代码省略.........