def residual_map_special_deltapsi_add_on( reflections,experiments,matches,hkllist, predicted,plot,eta_deg,deff ):

        detector = experiments[0].detector
        crystal = experiments[0].crystal
        unit_cell = crystal.get_unit_cell()
        pxlsz = detector[0].get_pixel_size()
        model_millers = reflections["miller_index"]
        dpsi = flex.double()
        for match in matches:

          obs_miller = hkllist[match["pred"]]
          model_index= model_millers.first_index(obs_miller)

          raw_delta_psi = reflections["delpsical.rad"][model_index]
          deltapsi_envelope = (unit_cell.d(obs_miller)/deff) + math.pi*eta_deg/180.
          normalized_delta_psi = raw_delta_psi/deltapsi_envelope

          dpsi.append( normalized_delta_psi )

        from matplotlib.colors import Normalize
        dnorm = Normalize()
        dnorm.autoscale(dpsi.as_numpy_array())

        CMAP = plot.get_cmap("bwr")
        for match,dcolor in zip(matches,dpsi):

          #print dcolor, dnorm(dcolor),  CMAP(dnorm(dcolor))
          #blue represents negative delta psi:  outside Ewald sphere; red, positive, inside Ewald sphere
          plot.plot([predicted[match["pred"]][1]/pxlsz[1]],[-predicted[match["pred"]][0]/pxlsz[0]],color=CMAP(dnorm(dcolor)),
          marker=".", markersize=5)

    def __init__(self, stretch='linear', exponent=5, vmid=None, vmin=None,
                 vmax=None, clip=False):
        '''
        Initalize an APLpyNormalize instance.

        Optional Keyword Arguments:

            *vmin*: [ None | float ]
                Minimum pixel value to use for the scaling.

            *vmax*: [ None | float ]
                Maximum pixel value to use for the scaling.

            *stretch*: [ 'linear' | 'log' | 'sqrt' | 'arcsinh' | 'power' ]
                The stretch function to use (default is 'linear').

            *vmid*: [ None | float ]
                Mid-pixel value used for the log and arcsinh stretches. If
                set to None, a default value is picked.

            *exponent*: [ float ]
                if self.stretch is set to 'power', this is the exponent to use.

            *clip*: [ True | False ]
                If clip is True and the given value falls outside the range,
                the returned value will be 0 or 1, whichever is closer.
        '''

        if vmax < vmin:
            raise Exception("vmax should be larger than vmin")

        # Call original initalization routine
        Normalize.__init__(self, vmin=vmin, vmax=vmax, clip=clip)

        # Save parameters
        self.stretch = stretch
        self.exponent = exponent

        if stretch == 'power' and np.equal(self.exponent, None):
            raise Exception("For stretch=='power', an exponent should be specified")

        if np.equal(vmid, None):
            if stretch == 'log':
                if vmin > 0:
                    self.midpoint = vmax / vmin
                else:
                    raise Exception("When using a log stretch, if vmin < 0, then vmid has to be specified")
            elif stretch == 'arcsinh':
                self.midpoint = -1./30.
            else:
                self.midpoint = None
        else:
            if stretch == 'log':
                if vmin < vmid:
                    raise Exception("When using a log stretch, vmin should be larger than vmid")
                self.midpoint = (vmax - vmid) / (vmin - vmid)
            elif stretch == 'arcsinh':
                self.midpoint = (vmid - vmin) / (vmax - vmin)
            else:
                self.midpoint = None

def model_to_pc2(model, x_start, y_start, resolution, width, height):
    """
    Creates a PointCloud2 by sampling a regular grid of points from the given model.
    """
    pc = PointCloud2()
    pc.header.stamp = rospy.get_rostime()
    pc.header.frame_id = 'map'
 
    xy_points = []
    for x in map_range(x_start, x_start + width, resolution):
        for y in map_range(y_start, y_start + height, resolution):
            xy_points.append([x, y])
    
    probs = model.score_samples(xy_points)
    
    # and normalise to range to make the visualisation prettier
    normaliser = Normalize()
    normaliser.autoscale(probs)
    probs = normaliser(probs)

    colour_map = plt.get_cmap('jet')    
    colours = colour_map(probs, bytes=True)

    cloud = []
    for i in range(len(probs)):
        cloud.append([xy_points[i][0], xy_points[i][1], 2*probs[i], pack_rgb(colours[i][0], colours[i][1], colours[i][2])])
 
    return create_cloud_xyzrgb(pc.header, cloud)

def visualize_temporal_activities(class_predictions, max_value=200, fps=1, title=None, legend=False):
    normalize = Normalize(vmin=1, vmax=max_value)
    normalize.clip=False
    cmap = plt.cm.Reds
    cmap.set_under('w')
    nb_instances = len(class_predictions)
    plt.figure(num=None, figsize=(18, 1), dpi=100)
    to_plot = class_predictions.astype(np.float32)
    to_plot[class_predictions==0.] = np.ma.masked
    plt.imshow(np.broadcast_to(to_plot, (20, nb_instances)), norm=normalize, interpolation='nearest', aspect='auto', cmap=cmap)
    if title:
        plt.title(title)
    ax = plt.gca()
    ax.get_yaxis().set_visible(False)

    if legend:
        index = np.arange(0,200)
        colors_index = np.unique(to_plot).astype(np.int64)
        if 0 in colors_index:
            colors_index = np.delete(colors_index, 0)
        patches = []
        for c in colors_index:
            patches.append(mpatches.Patch(color=cmap(normalize(c)), label=dataset.labels[c][1]))
        if patches:
            plt.legend(handles=patches, loc='upper center', bbox_to_anchor=(0.5, -.2), ncol=len(patches), fancybox=True, shadow=True)
    plt.show()

 def __init__(self):
     Normalize.__init__(self)
     self.stretch = "linear"
     self.bias = 0.5
     self.contrast = 0.5
     self.clip_lo = 5.0
     self.clip_hi = 95.0

 def __call__(self, value):
     if self.vmax <= self.vmin:
         self.vmax, self.vmin = self.vmin, self.vmax
         result = 1 - Normalize.__call__(self, value)
         self.vmax, self.vmin = self.vmin, self.vmax
     else:
         result = Normalize.__call__(self, value)
     return result

 def __init__(self, vmin=None, vmax=None, midpoint=0, clip=False):
     """
     Copied from SO answer: http://stackoverflow.com/questions/20144529/shifted-colorbar-matplotlib
     :param vmin:
     :param vmax:
     :param midpoint:
     :param clip:
     """
     self.midpoint = midpoint
     Normalize.__init__(self, vmin, vmax, clip)

def ImagePlot(image):
    if str(image.colorscale)=='n':
        remap = Normalize()
        remap.autoscale(image.data)
        ax.imshow(image.data, cmap='gray', norm=remap, origin='lower')
    elif str(image.colorscale) == 'yg':
        remap = LogNorm()
        remap.autoscale(image.data)
        ax.imshow(image.data, cmap='gray', norm=remap, origin='lower')
    elif str(image.colorscale) == 'ys':
        remap = LogNorm()
        remap.autoscale(image.data)
        ax.imshow(image.data, cmap='seismic', norm=remap, origin='lower')

def compare_temporal_activities(ground_truth, class_predictions, max_value=200, fps=1, title=None, legend=False, save_file='./img/activity_detection_sample_{}.png'):
    global count
    normalize = Normalize(vmin=1, vmax=max_value)
    normalize.clip=False
    cmap = plt.cm.Reds
    cmap.set_under('w')
    nb_instances = len(class_predictions)
    to_plot = np.zeros((20, nb_instances))
    to_plot[:10,:] = np.broadcast_to(ground_truth, (10, nb_instances))
    to_plot[10:,:] = np.broadcast_to(class_predictions, (10, nb_instances))
    to_plot = to_plot.astype(np.float32)
    to_plot[to_plot==0.] = np.ma.masked

    # Normalize the values and give them the largest distance possible between them
    unique_values = np.unique(to_plot).astype(np.int64)
    if 0 in unique_values:
        unique_values = np.delete(unique_values, 0)
    nb_different_values = len(unique_values)
    color_values = np.linspace(40, 190, nb_different_values)
    for i in range(nb_different_values):
        to_plot[to_plot == unique_values[i]] = color_values[i]

    plt.figure(num=None, figsize=(18, 1), dpi=100)
    plt.imshow(to_plot, norm=normalize, interpolation='nearest', aspect='auto', cmap=cmap)
    #plt.grid(True)
    plt.axhline(9, linestyle='-', color='k')
    plt.xlim([0,nb_instances])
    if title:
        plt.title(title)
    ax = plt.gca()
    #ax.get_yaxis().set_visible(False)
    ax.xaxis.grid(True, which='major')
    labels=['Ground\nTruth', 'Prediction']
    plt.yticks([5,15], labels, rotation="horizontal", size=13)
    plt.xlabel('Time (s)', horizontalalignment='left', fontsize=13)
    ax.xaxis.set_label_coords(0, -0.3)

    if legend:
        patches = []
        for c, l in zip(color_values, unique_values):
            patches.append(mpatches.Patch(color=cmap(normalize(c)), label=dataset.labels[l][1]))
        if patches:
            plt.legend(handles=patches, loc='upper center', bbox_to_anchor=(.5, -.2), ncol=len(patches), fancybox=True, shadow=True)
    #plt.show()
    plt.savefig(save_file.format(count), bbox_inches='tight')
    count += 1

def plot(d, sphere=False):
    """
    Plot directivity `d`.
    
    :param d: Directivity
    :type d: :class:`Directivity`
    
    :returns: Figure
    """
    
    #phi = np.linspace(-np.pi, +np.pi, 50)
    #theta = np.linspace(0.0, np.pi, 50)
    phi = np.linspace(0.0, +2.0*np.pi, 50)
    theta = np.linspace(0.0, np.pi, 50)
    THETA, PHI = np.meshgrid(theta, phi)
    
    # Directivity strength. Real-valued. Can be positive and negative.
    dr = d.using_spherical(THETA, PHI)
    
    if sphere:
        x, y, z = spherical_to_cartesian(1.0, THETA, PHI)
        
    else:
        x, y, z = spherical_to_cartesian( np.abs(dr), THETA, PHI )
    #R, THETA, PHI = cartesian_to_spherical(x, y, z)
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    #p = ax.plot_surface(x, y, z, cmap=plt.cm.jet, rstride=1, cstride=1, linewidth=0)


    norm = Normalize()
    norm.autoscale(dr)
    colors = cm.jet(norm(dr))
    m = cm.ScalarMappable(cmap=cm.jet, norm=norm)
    m.set_array(dr)
    p = ax.plot_surface(x, y, z, facecolors=colors, rstride=1, cstride=1, linewidth=0)
    plt.colorbar(m, ax=ax)
    
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$z$')
    return fig

    def __init__(self, vmin=None, vmax=None, midpoint=None,
                 clip=False, nlevs=9):
        """.. warning::
             MidpointNormalize is deprecated and will be removed in
             future versions.  Please use
             :class:`timutils.get_discrete_midpt_cmap_norm` instead

        returns a colormap and a matplotlib.colors.Normalize
        instance that implement a *continuous* colormap with an
        arbitrary midpoint.

        ARGS:
            vmin (real): the minimum value in the colormap
            vmax (real): the maximum value in the colormap
            midpoint (real): the midpoint to center on
            clip (boolean):
            nlevs (integer): number of levels to divide the colormap
                into.  Not currently functional.

        EXAMPLE:
            >>> import matplotlib.pyplot as plt
            >>> import numpy as np
            >>> from timutils.midpt_norm import MidpointNormalize
            >>> plt.close('all')
            >>> data = np.random.randint(-120, 20, [124, 124])
            >>> fix, ax = plt.subplots(1, 2)
            >>> mycmap = plt.get_cmap('Blues')
            >>> mynorm = MidpointNormalize(vmin=-120, vmax=20, midpoint=0.0)
            >>> cm = ax[0].pcolormesh(data, norm=mynorm, cmap=mycmap)
            >>> plt.colorbar(cm, cax=ax[1], norm=mynorm, cmap=mycmap)
            >>> plt.show()

        adapted by Timothy W. Hilton from `code posted by Joe Kington
        <http://stackoverflow.com/questions/20144529/shifted-colorbar-matplotlib>`_
        accessed 19 January 2015
        """
        warnings.warn(('MidpointNormalize is (1) deprecated and (2)'
                       'buggy and will be' 'removed in future versions. '
                       'Please use get_discrete_midpt_cmap_norm instead'))
        self.midpoint = midpoint
        self.nlevs = nlevs
        Normalize.__init__(self, vmin, vmax, clip)

    def auto_scale_cross_plot(self, event):
        
        norm = Normalize()
        
        for hl in self.h_cross_slice_plot.get_lines(): 
            d = hl.get_ydata()
            norm.autoscale(d)
            hl.set_ydata(norm(d))
          
        for vl in self.v_cross_slice_plot.get_lines(): 
            d = vl.get_ydata()
            norm.autoscale(d)
            vl.set_ydata(norm(d))
        
        
        self.v_cross_slice_plot.relim()
        self.h_cross_slice_plot.relim()
        self.v_cross_slice_plot.autoscale_view(True,True,True)
        self.h_cross_slice_plot.autoscale_view(True,True,True)

        self.cross_slice_canvas.draw()

  def plot(self,dano_summation):
    from matplotlib import pyplot as plt

    if self.params.use_weights:
      wt = 1./(self.diffs.sigmas()*self.diffs.sigmas())
      order = flex.sort_permutation(wt)
      wt = wt.select(order)
      df = self.diffs.data().select(order)
      dano = dano_summation.select(self.sel0).select(order)
      from matplotlib.colors import Normalize
      dnorm = Normalize()
      dnorm.autoscale(wt.as_numpy_array())
      CMAP = plt.get_cmap("rainbow")
      for ij in xrange(len(self.diffs.data())):
        #blue represents zero weight:  red, large weight
        plt.plot([df[ij]],[dano[ij]],color=CMAP(dnorm(wt[ij])),marker=".", markersize=4)

    else:
      plt.plot(self.diffs.data(),dano_summation.select(self.sel0),"r,")
    plt.axes().set_aspect("equal")
    plt.axes().set_xlabel("Observed Dano")
    plt.axes().set_ylabel("Model Dano")
    plt.show()

def plot2D(X, filename=None, last_column_color=False):
    x1 = X[:, 0]
    x2 = X[:, 1]
    m = X.shape[0]
    if last_column_color:
        c = X[:, -1]
        c_map = get_cmap('jet')
        c_norm = Normalize()
        c_norm.autoscale(c)
        scalar_map = ScalarMappable(norm=c_norm, cmap=c_map)
        color_val = scalar_map.to_rgba(c)
    else:
        color_val = 'b' * m
    fig = figure()
    ax = fig.add_subplot(111)
    for i in range(m):
        ax.plot(x1[i], x2[i], 'o', color=color_val[i])
    if filename is None:
        fig.show()
    else:
        fig.savefig(filename + ".png")
    fig.clf()
    close()

 def __init__(self,vmin=None,vmax=None,clip=False):
     Normalize.__init__(self,vmin,vmax,clip)

 def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
     self.midpoint = midpoint
     # m=max(abs(vmin),abs(vmax))
     Normalize.__init__(self, vmin, vmax, clip)

 def __init__(self,linthresh,vmin=None,vmax=None,clip=False):
   Normalize.__init__(self,vmin,vmax,clip)
   self.linthresh=linthresh
   self.vmin, self.vmax = vmin, vmax

 def __call__(self, epoch, **kargs):
     if isinstance(epoch, np.ndarray) and len(epoch.shape) == 0:
         epoch = epoch.item()
     return Normalize.__call__(self, self.mktime(epoch), **kargs)

 def __init__(self, epoch0, epoch1):
     Normalize.__init__(self, self.mktime(epoch0), self.mktime(epoch1))

  def plot_one_model(self,nrow,out):
    fig = plt.subplot(self.gs[nrow*self.ncols])
    two_thetas = self.reduction.get_two_theta_deg()
    degrees = self.reduction.get_delta_psi_deg()

    if self.color_encoding=="conventional":
          positive = (self.reduction.i_sigi>=0.)
          fig.plot(two_thetas.select(positive), degrees.select(positive), "bo")
          fig.plot(two_thetas.select(~positive), degrees.select(~positive), "r+")
    elif self.color_encoding=="I/sigma":
          positive = (self.reduction.i_sigi>=0.)
          tt_selected = two_thetas.select(positive)
          dp_selected = degrees.select(positive)
          i_sigi_select = self.reduction.i_sigi.select(positive)
          order = flex.sort_permutation(i_sigi_select)
          tt_selected = tt_selected.select(order)
          dp_selected = dp_selected.select(order)
          i_sigi_selected = i_sigi_select.select(order)
          from matplotlib.colors import Normalize
          dnorm = Normalize()
          dcolors = i_sigi_selected.as_numpy_array()
          dnorm.autoscale(dcolors)
          N = len(dcolors)
          CMAP = plt.get_cmap("rainbow")
          if self.refined.get("partiality_array",None) is None:
            for n in xrange(N):
              fig.plot([tt_selected[n]],[dp_selected[n]],
              color=CMAP(dnorm(dcolors[n])),marker=".", markersize=10)
          else:
            partials = self.refined.get("partiality_array")
            partials_select = partials.select(positive)
            partials_selected = partials_select.select(order)
            assert len(partials)==len(positive)
            for n in xrange(N):
              fig.plot([tt_selected[n]],[dp_selected[n]],
              color=CMAP(dnorm(dcolors[n])),marker=".", markersize=20*partials_selected[n])
              # change the markersize to indicate partiality.
          negative = (self.reduction.i_sigi<0.)
          fig.plot(two_thetas.select(negative), degrees.select(negative), "r+", linewidth=1)
    else:
          strong = (self.reduction.i_sigi>=10.)
          positive = ((~strong) & (self.reduction.i_sigi>=0.))
          negative = (self.reduction.i_sigi<0.)
          assert (strong.count(True)+positive.count(True)+negative.count(True) ==
                  len(self.reduction.i_sigi))
          fig.plot(two_thetas.select(positive), degrees.select(positive), "bo")
          fig.plot(two_thetas.select(strong), degrees.select(strong), marker='.',linestyle='None',
           markerfacecolor='#00ee00', markersize=10)
          fig.plot(two_thetas.select(negative), degrees.select(negative), "r+")

    # indicate the imposed resolution filter
    wavelength = self.reduction.experiment.beam.get_wavelength()
    imposed_res_filter = self.reduction.get_imposed_res_filter(out)
    resolution_markers = [
      a for a in [imposed_res_filter,self.reduction.measurements.d_min()] if a is not None]
    for RM in resolution_markers:
          two_th = (180./math.pi)*2.*math.asin(wavelength/(2.*RM))
          plt.plot([two_th, two_th],[self.AD1TF7B_MAXDP*-0.8,self.AD1TF7B_MAXDP*0.8],'k-')
          plt.text(two_th,self.AD1TF7B_MAXDP*-0.9,"%4.2f"%RM)

    #indicate the linefit
    mean = flex.mean(degrees)
    minplot = flex.min(two_thetas)
    plt.plot([0,minplot],[mean,mean],"k-")
    LR = flex.linear_regression(two_thetas, degrees)
    model_y = LR.slope()*two_thetas + LR.y_intercept()
    plt.plot(two_thetas, model_y, "k-")

    #Now let's take care of the red and green lines.
    half_mosaic_rotation_deg = self.refined["half_mosaic_rotation_deg"]
    mosaic_domain_size_ang = self.refined["mosaic_domain_size_ang"]
    red_curve_domain_size_ang = self.refined.get("red_curve_domain_size_ang",mosaic_domain_size_ang)
    a_step = self.AD1TF7B_MAX2T / 50.
    a_range = flex.double([a_step*x for x in xrange(1,50)]) # domain two-theta array
    #Bragg law [d=L/2sinTH]
    d_spacing = (wavelength/(2.*flex.sin(math.pi*a_range/360.)))
    # convert two_theta to a delta-psi.  Formula for Deffective [Dpsi=d/2Deff]
    inner_phi_deg = flex.asin((d_spacing / (2.*red_curve_domain_size_ang)) )*(180./math.pi)
    outer_phi_deg = flex.asin((d_spacing / (2.*mosaic_domain_size_ang)) + \
      half_mosaic_rotation_deg*math.pi/180. )*(180./math.pi)
    plt.title("ML: mosaicity FW=%4.2f deg, Dsize=%5.0fA on %d spots\n%s"%(
          2.*half_mosaic_rotation_deg, mosaic_domain_size_ang, len(two_thetas),
          os.path.basename(self.reduction.filename)))
    plt.plot(a_range, inner_phi_deg, "r-")
    plt.plot(a_range,-inner_phi_deg, "r-")
    plt.plot(a_range, outer_phi_deg, "g-")
    plt.plot(a_range, -outer_phi_deg, "g-")
    plt.xlim([0,self.AD1TF7B_MAX2T])
    plt.ylim([-self.AD1TF7B_MAXDP,self.AD1TF7B_MAXDP])

    #second plot shows histogram
    fig = plt.subplot(self.gs[1+nrow*self.ncols])
    plt.xlim([-self.AD1TF7B_MAXDP,self.AD1TF7B_MAXDP])
    nbins = 50
    n,bins,patches = plt.hist(dp_selected, nbins,
           range=(-self.AD1TF7B_MAXDP,self.AD1TF7B_MAXDP),
           weights=self.reduction.i_sigi.select(positive),
           normed=0, facecolor="orange", alpha=0.75)
    #ersatz determine the median i_sigi point:
    isi_positive = self.reduction.i_sigi.select(positive)
#.........这里部分代码省略.........