def __call__(self, transform_xy, x1, y1, x2, y2):
        """
        get extreme values.

        x1, y1, x2, y2 in image coordinates (0-based)
        nx, ny : number of divisions in each axis
        """
        x_, y_ = np.linspace(x1, x2, self.nx), np.linspace(y1, y2, self.ny)
        x, y = np.meshgrid(x_, y_)
        lon, lat = transform_xy(np.ravel(x), np.ravel(y))

        # iron out jumps, but algorithm should be improved.
        # This is just naive way of doing and my fail for some cases.
        # Consider replacing this with numpy.unwrap
        # We are ignoring invalid warnings. They are triggered when
        # comparing arrays with NaNs using > We are already handling
        # that correctly using np.nanmin and np.nanmax
        with np.errstate(invalid='ignore'):
            if self.lon_cycle is not None:
                lon0 = np.nanmin(lon)
                lon -= 360. * ((lon - lon0) > 180.)
            if self.lat_cycle is not None:
                lat0 = np.nanmin(lat)
                lat -= 360. * ((lat - lat0) > 180.)

        lon_min, lon_max = np.nanmin(lon), np.nanmax(lon)
        lat_min, lat_max = np.nanmin(lat), np.nanmax(lat)

        lon_min, lon_max, lat_min, lat_max = \
                 self._adjust_extremes(lon_min, lon_max, lat_min, lat_max)

        return lon_min, lon_max, lat_min, lat_max

    def _set_minmax(self):
        data = self._get_fast_data()
        try:
            self.maxval = numpy.nanmax(data)
            self.minval = numpy.nanmin(data)
        except Exception:
            self.maxval = 0
            self.minval = 0

        # TODO: see if there is a faster way to ignore infinity
        try:
            if numpy.isfinite(self.maxval):
                self.maxval_noinf = self.maxval
            else:
                self.maxval_noinf = numpy.nanmax(data[numpy.isfinite(data)])
        except:
            self.maxval_noinf = self.maxval

        try:
            if numpy.isfinite(self.minval):
                self.minval_noinf = self.minval
            else:
                self.minval_noinf = numpy.nanmin(data[numpy.isfinite(data)])
        except:
            self.minval_noinf = self.minval

def summary():
    # read sonde data
    for sites in [[0],[1],[2]]:
        slist,snames=read_diff_events(sites=sites)
        ecount = [len(s.einds) for s in slist]
        mintp = [np.nanmin(s.tp) for s in slist]
        meantp = [np.nanmean(s.tp) for s in slist]
        maxtp = [np.nanmax(s.tp) for s in slist]
        
        head="%9s"%slist[0].name
        ecount = "events   "
        meantp = "mean tph "
        minmax = "tph bound"
        for sonde, sname in zip(slist,snames):
            
            head=head+'| %16s'%sname
            ecount=ecount+'| %16d'%len(sonde.einds)
            meantp=meantp+'| %16.2f'%np.nanmean(sonde.tp)
            minmax=minmax+'| %7.2f,%7.2f '%(np.nanmin(sonde.tp),np.nanmax(sonde.tp))
            
        print("")
        print(head)
        print(ecount)
        print(meantp)
        print(minmax)

    def _get_Tp_limits(self):
        """Get the limits for the graphs in temperature and pressure, based on 
        SI units: [Tmin, Tmax, pmin, pmax]"""
        T_lo,T_hi,P_lo,P_hi = self.limits
        Ts_lo,Ts_hi = self._get_sat_bounds(CoolProp.iT)
        Ps_lo,Ps_hi = self._get_sat_bounds(CoolProp.iP)

        if T_lo is None:            T_lo  = 0.0
        elif T_lo < self.ID_FACTOR: T_lo *= Ts_lo
        if T_hi is None:            T_hi  = 1e6
        elif T_hi < self.ID_FACTOR: T_hi *= Ts_hi
        if P_lo is None:            P_lo  = 0.0
        elif P_lo < self.ID_FACTOR: P_lo *= Ps_lo
        if P_hi is None:            P_hi  = 1e10
        elif P_hi < self.ID_FACTOR: P_hi *= Ps_hi

        try: T_lo = np.nanmax([T_lo, self._state.trivial_keyed_output(CoolProp.iT_min)])
        except: pass
        try: T_hi = np.nanmin([T_hi, self._state.trivial_keyed_output(CoolProp.iT_max)])
        except: pass
        try: P_lo = np.nanmax([P_lo, self._state.trivial_keyed_output(CoolProp.iP_min)])
        except: pass
        try: P_hi = np.nanmin([P_hi, self._state.trivial_keyed_output(CoolProp.iP_max)])
        except: pass

        return [T_lo,T_hi,P_lo,P_hi]

def calc_norm_summary_tables(accuracy_tbl, time_tbl):
    """
    Calculate normalized performance/ranking summary, as numpy
    matrices as usual for convenience, and matrices of additional
    statistics (min, max, percentiles, etc.)

    Here normalized means relative to the best which gets a 1, all
    others get the ratio resulting from dividing by the performance of
    the best.
    """
    # Min across all minimizers, i.e. for each fit problem what is the lowest chi-squared and the lowest time
    min_sum_err_sq = np.nanmin(accuracy_tbl, 1)
    min_runtime = np.nanmin(time_tbl, 1)

    # create normalised tables
    norm_acc_rankings = accuracy_tbl / min_sum_err_sq[:, None]
    norm_runtimes = time_tbl / min_runtime[:, None]

    summary_cells_acc = np.array([np.nanmin(norm_acc_rankings, 0),
                                  np.nanmax(norm_acc_rankings, 0),
                                  stats.nanmean(norm_acc_rankings, 0),
                                  stats.nanmedian(norm_acc_rankings, 0)
                                  ])

    summary_cells_runtime = np.array([np.nanmin(norm_runtimes, 0),
                                      np.nanmax(norm_runtimes, 0),
                                      stats.nanmean(norm_runtimes, 0),
                                      stats.nanmedian(norm_runtimes, 0)
                                      ])

    return norm_acc_rankings, norm_runtimes, summary_cells_acc, summary_cells_runtime

def classify(request):
    C = json.loads(request.POST["C"])
    try:
        features, labels = get_multi_features(request)
    except ValueError as e:
        return HttpResponse(json.dumps({"status": e.message}))
    try:
        kernel = get_kernel(request, features)
    except ValueError as e:
        return HttpResponse(json.dumps({"status": e.message}))
    
    learn = "No"  
    values=[]

    try:
        domain = json.loads(request.POST['axis_domain'])
        x, y, z = svm.classify_svm(sg.GMNPSVM, features, labels, kernel, domain, learn, values, C, False)
    except Exception as e:
        return HttpResponse(json.dumps({"status": repr(e)}))

#    z = z + np.random.rand(*z.shape) * 0.01
	
    z_max = np.nanmax(z)
    z_min = np.nanmin(z)
    z_delta = 0.1*(np.nanmax(z)-np.nanmin(z))
    data = {"status": "ok",
            "domain": [z_min-z_delta, z_max+z_delta],
            "max": z_max+z_delta,
            "min": z_min-z_delta,
            "z": z.tolist()}

    return HttpResponse(json.dumps(data))

	def show(self,**kwargs):
		display = kwargs.get('display', True)
		show_layers = kwargs.get('show_layers',self.layers)
		try:
			show_layers=sorted(show_layers)
		except TypeError:
			show_layers=[show_layers]
		extent=kwargs.get('extent', 
						max_axis(*tuple(_image.axis for _image in self.image_sorted[self.layers[0]])))
		vmin=kwargs.get('vmin')
		vmax=kwargs.get('vmax')
		fig = plt.figure(figsize=(8, 8*abs((extent[3]-extent[2])*1./(extent[1]-extent[0]))))
		for layer in show_layers:
			for image in self.image_sorted[layer]:
				if layer==show_layers[0] and image==self.image_sorted[layer][0]:
					if not vmin:
						kwargs['vmin']=np.nanmin(image.image)
						vmin=np.nanmin(image.image)
					if not vmax:
						kwargs['vmax']=np.nanmax(image.image)
						vmax=np.nanmax(image.image)
					image.show(hold=True,**kwargs)
				else:
					image.show(hold=True,vmin=vmin,vmax=vmax,scalebar='off',colorbar='off')
		plt.xlim(extent[:2])
		plt.ylim(extent[-2:])
		if display:
			plt.show()
		else:
			return fig

def bin_fit(x, y, buckets=3):
     
    assert buckets in [3,25]

    xstd=np.nanstd(x)
    
    if buckets==3:
        binlimits=[np.nanmin(x), -xstd/2.0,xstd/2.0 , np.nanmax(x)]
    elif buckets==25:
    
        steps=xstd/4.0
        binlimits=np.arange(-xstd*3.0, xstd*3.0, steps)
    
        binlimits=[np.nanmin(x)]+list(binlimits)+[np.nanmax(x)]
    
    fit_y=[]
    err_y=[]
    x_values_to_plot=[]
    for binidx in range(len(binlimits))[1:]:
        lower_bin_x=binlimits[binidx-1]
        upper_bin_x=binlimits[binidx]

        x_values_to_plot.append(np.mean([lower_bin_x, upper_bin_x]))

        y_in_bin=[y[idx] for idx in range(len(y)) if x[idx]>=lower_bin_x and x[idx]<upper_bin_x]

        fit_y.append(np.nanmedian(y_in_bin))
        err_y.append(np.nanstd(y_in_bin))

    ## no zeros
    

    return (binlimits, x_values_to_plot, fit_y, err_y)

 def test_threshold_filter_nan(self):
     src = self.make_src(nan=True)
     self.e.add_source(src)
     threshold = Threshold()
     self.e.add_filter(threshold)
     self.assertEqual(np.nanmin(src.scalar_data), np.nanmin(threshold.outputs[0].point_data.scalars.to_array()))
     self.assertEqual(np.nanmax(src.scalar_data), np.nanmax(threshold.outputs[0].point_data.scalars.to_array()))

def TestPlot(fig=None):
    A = numpy.array([1,2,3,4,2,5,8,3,2,3,5,6])
    B = numpy.array([8,7,3,6,4,numpy.nan,9,3,7,numpy.nan,2,4])
    C = numpy.array([6,3,4,7,2,1,1,7,8,4,3,2])
    D = numpy.array([5,2,4,5,3,8,2,5,3,5,6,8])
    
    # A work around to get the histograms overplotted with each other to overlap correctly;
    histrangelist = [(numpy.nanmin(A),numpy.nanmax(A)),(numpy.nanmin(B),numpy.nanmax(B)),
                (numpy.nanmin(C),numpy.nanmax(C)),(numpy.nanmin(D),numpy.nanmax(D))]
    
    data = numpy.array([A,B,C,D])
    labels = ['A','3','C','D']

    fig = GridPlot(data,labels=labels, no_tick_labels=True, color='black', 
                    hist=True, histbins=3, histloc='tl', histrangelist=histrangelist, fig=None) 
    
    # Data of note to plot in different color
    A2 = numpy.array([1,2,3,4])
    B2 = numpy.array([8,7,3,6])
    C2 = numpy.array([6,3,4,7])
    D2 = numpy.array([5,2,4,5])
    data2 = numpy.array([A2,B2,C2,D2])
    
    fig = GridPlot(data2,labels=labels, no_tick_labels=True, color='red', 
                hist=True, histbins=3, histloc='tr', histrangelist=histrangelist, fig=fig) 
    
    return fig

    def plot_result(self, result):
        """
        It plots the resulting Q and q when atype is set to 'tsl' or 'asl'

         :param result:
           Event Sync result from compute()
        :type result: dict

        :returns: plt.figure
               -- figure plot
        """

        ' Raise error if parameters are not in the correct type '
        if not(isinstance(result, dict)) : raise TypeError("Requires result to be a dictionary")

        ' Raise error if not the good dictionary '
        if not 'Q' in result : raise ValueError("Requires dictionary to be the output of compute() method")
        if not 'q' in result : raise ValueError("Requires dictionary to be the output of compute() method")

        x=np.arange(0, result['Q'].size, 1)

        figure, axarr = plt.subplots(2, sharex=True)
        axarr[0].set_title('Synchrony and time delay pattern')
        axarr[0].set_xlabel('Samples')
        axarr[1].set_xlabel('Samples')
        axarr[0].set_ylim(0,np.nanmax(result['Q']))
        axarr[0].plot(x, result['Q'], label="Synchrony (Qn)")
        axarr[1].set_ylim(np.nanmin(result['q']),np.nanmax(result['q']))
        axarr[1].plot(x, result['q'], label="Time delay pattern (qn)")
        axarr[0].legend(loc='best')
        axarr[1].legend(loc='best')

        return figure

def plot_all_time_series(config_list, output_dir):
    """Plot column charts of the raw total time/energy spent in each profiler category.

    Keyword arguments:
    config_list -- [(config, result of process_config_dir(...))]
    output_dir -- where to write plots to
    """
    time_series_out_dir = path.join(output_dir, 'time_series')
    os.makedirs(time_series_out_dir)

    max_end_times = []
    max_power_values = []
    for (c, cd) in config_list:
        for (t, td) in cd:
            trial_max_end_times = map(np.nanmax, filter(lambda x: len(x) > 0, [te for (p, ts, te, es, ee) in td]))
            max_end_times.append(np.nanmax(trial_max_end_times))
            for (p, ts, te, es, ee) in td:
                # We only care about the energy profiler (others aren't reliable for instant power anyway)
                if p == ENERGY_PROFILER_NAME and len(te) > 0:
                    max_power_values.append(np.nanmax(hb_energy_times_to_power(es, ee, ts, te)))
    max_time = np.nanmax(max_end_times)
    max_power = np.nanmax(np.array(max_power_values)) * 1.2  # leave a little space at the top

    for (config, config_data) in config_list:
        [plot_trial_time_series(config, trial, trial_data, max_time, max_power, time_series_out_dir)
            for (trial, trial_data) in config_data]

def plot_nontarget_betas_n_back(t_vols_n_back_beta_1, b_vols_smooth_n_back, in_brain_mask, brain_structure, nice_cmap, n_back):

  beta_index = 1

  b_vols_smooth_n_back[~in_brain_mask] = np.nan
  t_vols_n_back_beta_1[~in_brain_mask] = np.nan
  min_val = np.nanmin(b_vols_smooth_n_back[...,(40,50,60),beta_index])
  max_val = np.nanmax(b_vols_smooth_n_back[...,(40,50,60),beta_index])

  plt.figure()

  for map_index, depth in (((3,2,1), 40),((3,2,3), 50),((3,2,5), 60)):
    plt.subplot(*map_index)
    plt.title("z=%d,%s" % (depth, n_back + "-back nontarget,beta values"))
    plt.imshow(brain_structure[...,depth], alpha=0.5)
    plt.imshow(b_vols_smooth_n_back[...,depth,beta_index], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
    plt.colorbar()
    plt.tight_layout()

  t_min_val = np.nanmin(t_vols_n_back_beta_1[...,(40,50,60)])
  t_max_val = np.nanmax(t_vols_n_back_beta_1[...,(40,50,60)])

  for map_index, depth in (((3,2,2), 40),((3,2,4), 50),((3,2,6), 60)):
    plt.subplot(*map_index)
    plt.title("z=%d,%s" % (depth, n_back + "-back nontarget,t values"))
    plt.imshow(brain_structure[...,depth], alpha=0.5)
    plt.imshow(t_vols_n_back_beta_1[...,depth], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
    plt.colorbar()
    plt.tight_layout()

  plt.savefig(os.path.join(output_filename, "sub011_nontarget_betas_%s_back.png" % (n_back)), format='png', dpi=500)  

def plot_noise_regressor_betas(b_vols_smooth, t_vols_beta_6_to_9, brain_structure, in_brain_mask, nice_cmap):

  plt.figure()

  min_val = np.nanmin(b_vols_smooth[...,40,(6,7,9)])
  max_val = np.nanmax(b_vols_smooth[...,40,(6,7,9)])

  plt.subplot(3,2,1)
  plt.title("z=%d,%s" % (40, "linear drift,betas"))
  b_vols_smooth[~in_brain_mask] = np.nan
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(b_vols_smooth[...,40,6], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,3)
  plt.title("z=%d,%s" % (40, "quadratic drift,betas"))
  b_vols_smooth[~in_brain_mask] = np.nan
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(b_vols_smooth[...,40,7], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,5)
  plt.title("z=%d,%s" % (40, "second PC,betas"))
  b_vols_smooth[~in_brain_mask] = np.nan
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(b_vols_smooth[...,40,9], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
  plt.colorbar()
  plt.tight_layout()

  t_vols_beta_6_to_9[0][~in_brain_mask] = np.nan
  t_vols_beta_6_to_9[1][~in_brain_mask] = np.nan
  t_vols_beta_6_to_9[3][~in_brain_mask] = np.nan

  t_min_val = np.nanmin([t_vols_beta_6_to_9[0][...,40], t_vols_beta_6_to_9[1][...,40], t_vols_beta_6_to_9[3][...,40]])
  t_max_val = np.nanmax([t_vols_beta_6_to_9[0][...,40], t_vols_beta_6_to_9[1][...,40], t_vols_beta_6_to_9[3][...,40]])

  plt.subplot(3,2,2)
  plt.title("z=%d,%s" % (40, "linear drift,t values"))
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(t_vols_beta_6_to_9[0][...,40], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,4)
  plt.title("z=%d,%s" % (40, "quadratic drift,t values"))
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(t_vols_beta_6_to_9[1][...,40], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,6)
  plt.title("z=%d,%s" % (40, "second PC,t values"))
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(t_vols_beta_6_to_9[3][...,40], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.savefig(os.path.join(output_filename, "sub001_noise_regressors_betas_map.png"), format='png', dpi=500)  

    def acquire_data(self, var_name=None, slice_=()):
        if var_name in self._variables:
            vars = [var_name]
        else:
            vars = self._variables

        if not isinstance(slice_, tuple): slice_ = (slice_,)

        for vn in vars:
            var = self._data_array[vn]

            ndims = len(var.shape)
            # Ensure the slice_ is the appropriate length
            if len(slice_) < ndims:
                slice_ += (slice(None),) * (ndims-len(slice_))

            arri = ArrayIterator(var, self._block_size)[slice_]
            for d in arri:
                if d.dtype.char is "S":
                    # Obviously, we can't get the range of values for a string data type!
                    rng = None
                elif isinstance(d, numpy.ma.masked_array):
                    # TODO: This is a temporary fix because numpy 'nanmin' and 'nanmax'
                    # are currently broken for masked_arrays:
                    # http://mail.scipy.org/pipermail/numpy-discussion/2011-July/057806.html
                    dc = d.compressed()
                    if dc.size == 0:
                        rng = None
                    else:
                        rng = (numpy.nanmin(dc), numpy.nanmax(dc))
                else:
                    rng = (numpy.nanmin(d), numpy.nanmax(d))
                yield vn, arri.curr_slice, rng, d

        return

def plot_richness_scatter(gals, name, full_set):
    log_counts_a, scatter_a = richness_scatter(gals[gals['ssfr'] < -11.0], full_set)
    log_counts_p, scatter_p = richness_scatter(gals[gals['pred'] < -11.0], full_set)
    #fig1 = plt.figure(figsize=(12,7))
    #frame1=fig1.add_axes((.1,.3,.8,.6))
    #plt.subplot(121)
    plt.plot(log_counts_a, scatter_a, 'o', label='input', color='k', markersize=7)
    plt.plot(log_counts_p, scatter_p, 'o', label='predicted', color=red_col, markersize=7)
    #plt.title('Scatter in richness ' + name)
    plt.xlabel('Log Number of red satellites')
    plt.xlabel('$<log N_{red sat}>$')
    plt.xlim(-.1,2.6)
    plt.ylim(0, np.max([np.nanmax(scatter_a),np.nanmax(scatter_p)]) +.1)
    plt.ylabel('Scatter in $M_{halo}$')
    plt.legend(loc='best')

    #plt.subplot(122)
    #frame2=fig1.add_axes((.1,.1,.8,.2))
    #series_a = pd.Series(scatter_a, index=counts_a)
    #series_p = pd.Series(scatter_p, index=counts_p)
    # scat_diff = (series_a - series_p)/series_a
    #scat_ratio = series_p/series_a
    #plt.plot(scat_diff.index, scat_diff.values, 'ob')
    #plt.plot(scat_ratio.index, scat_ratio.values, 'ob')
    #plt.title("Scatter ratios in richness for actual vs predicted")
    #plt.axhline(0)
    #plt.ylabel('Error')
    #plt.xlabel('Number of red satellites')
    return

 def preprocess_crowdlabels(self, crowdlabels):
     # Initialise all objects relating to the crowd labels.
     C = {}
     crowdlabels[np.isnan(crowdlabels)] = -1
     if self.discretedecisions:
         crowdlabels = np.round(crowdlabels).astype(int)
     if self.table_format_flag:# crowd labels as a full KxN table? If false, use diags sparse 3-column list, where 1st
         logging.error("Can't use table format with preference pairs at the moment.")
         return
     
     if self.K < int(np.nanmax(crowdlabels[:,0]))+1:
         self.K = int(np.nanmax(crowdlabels[:,0]))+1 # add one because indexes start from 0
         
     for l in range(self.nscores):
         lIdxs = np.argwhere(crowdlabels[:, 3] == l)[:,0]
         data = np.ones((len(lIdxs), 1)).reshape(-1)
         rows = np.array(crowdlabels[lIdxs, 1]).reshape(-1) * self.N + crowdlabels[lIdxs, 2]
         cols = np.array(crowdlabels[lIdxs, 0]).reshape(-1)
         
         Cl = csr_matrix(coo_matrix((data,(rows,cols)), shape=(self.N**2, self.K)))
         C[l] = Cl            
         
     # Set and reset object properties for the new dataset
     self.C = C
     self.lnpCT = np.zeros((self.N, self.nclasses))
     self.conf_mat_ind = []
     # pre-compute the indices into the pi arrays
     # repeat for test labels only
     
     self.Ctest = {}
     for l in range(self.nscores):
         self.Ctest[l] = C[l][self.testidxs, :]
     # Reset the pre-calculated data for the training set in case goldlabels has changed
     self.alpha_tr = []

def draw_hmap_old(hmap, yvals, fname=None):
    """
    Plot a matrix as a heat map and write an image file.
    :param hmap: Heat map matrix.
    :param yvals: Heat map Y labels (e.g. amino acid names).
    :param fname: Destination image file.
    """
    if np.nanmax(hmap) > abs(np.nanmin(hmap)):
        vmax = np.nanmax(hmap)
        vmin = -np.nanmax(hmap)
    else:
        vmax = abs(np.nanmin(hmap))
        vmin = np.nanmin(hmap)
    fig = plt.figure()
    plt.figure(figsize=(20,10))
    plt.imshow(hmap, cmap='RdBu', interpolation = 'nearest',aspect='auto',vmin = vmin ,vmax = vmax )
    plt.xlim(0, hmap.shape[1])
    plt.ylim(0, hmap.shape[0])
    ax = plt.gca()
    fig.set_facecolor('white')
    ax.set_xlim((-0.5, hmap.shape[1] -0.5))
    ax.set_ylim((-0.5, hmap.shape[0] -0.5))
    ax.set_yticks([x for x in xrange(0, hmap.shape[0])])
    ax.set_yticklabels(yvals)
    ax.set_xticks(range(0,76,5))
    ax.set_xticklabels(range(2,76,5)+['STOP'])
    ax.set_ylabel('Residue')
    ax.set_xlabel('Ub Sequence Position')
    cb = plt.colorbar()
    cb.set_clim(vmin=vmin, vmax=vmax)
    cb.set_label('Relative Fitness')
    if fname is not None:
        plt.savefig(fname, bbox_inches='tight')
    return fig

    def _axes_domain(self, nx=None, ny=None, background_patch=None):
        """Returns x_range, y_range"""
        DEBUG = False

        transform = self._crs_transform()

        ax_transform = self.axes.transAxes
        desired_trans = ax_transform - transform

        nx = nx or 30
        ny = ny or 30
        x = np.linspace(1e-9, 1 - 1e-9, nx)
        y = np.linspace(1e-9, 1 - 1e-9, ny)
        x, y = np.meshgrid(x, y)

        coords = np.concatenate([x.flatten()[:, None], y.flatten()[:, None]], 1)

        in_data = desired_trans.transform(coords)

        ax_to_bkg_patch = self.axes.transAxes - background_patch.get_transform()

        ok = np.zeros(in_data.shape[:-1], dtype=np.bool)
        # XXX Vectorise contains_point
        for i, val in enumerate(in_data):
            # convert the coordinates of the data to the background
            # patches coordinates
            background_coord = ax_to_bkg_patch.transform(coords[i : i + 1, :])
            bkg_patch_contains = background_patch.get_path().contains_point
            if bkg_patch_contains(background_coord[0, :]):
                color = "r"
                ok[i] = True
            else:
                color = "b"

            if DEBUG:
                import matplotlib.pyplot as plt

                plt.plot(coords[i, 0], coords[i, 1], "o" + color, clip_on=False, transform=ax_transform)
        #                plt.text(coords[i, 0], coords[i, 1], str(val), clip_on=False,
        #                         transform=ax_transform, rotation=23,
        #                         horizontalalignment='right')

        inside = in_data[ok, :]
        x_range = np.nanmin(inside[:, 0]), np.nanmax(inside[:, 0])
        y_range = np.nanmin(inside[:, 1]), np.nanmax(inside[:, 1])

        # XXX Cartopy specific thing. Perhaps make this bit a specialisation
        # in a subclass...
        crs = self.crs
        if isinstance(crs, Projection):
            x_range = np.clip(x_range, *crs.x_limits)
            y_range = np.clip(y_range, *crs.y_limits)

            # if the limit is >90 of the full x limit, then just use the full
            # x limit (this makes circular handling better)
            prct = np.abs(np.diff(x_range) / np.diff(crs.x_limits))
            if prct > 0.9:
                x_range = crs.x_limits

        return x_range, y_range

 def test_derivative(self):
     log= logging.getLogger( "test_J0023.derivative_test")
     testp = tdu.get_derivative_params(self.modelJ0023)
     delay = self.modelJ0023.delay(self.toasJ0023)
     for p in testp.keys():
         log.debug( "Runing derivative for %s", 'd_delay_d_'+p)
         if p in ['EPS2', 'EPS1']:
             testp[p] = 10
         ndf = self.modelJ0023.d_phase_d_param_num(self.toasJ0023, p, testp[p])
         adf = self.modelJ0023.d_phase_d_param(self.toasJ0023, delay, p)
         diff = adf - ndf
         if not np.all(diff.value) == 0.0:
             mean_der = (adf+ndf)/2.0
             relative_diff = np.abs(diff)/np.abs(mean_der)
             #print "Diff Max is :", np.abs(diff).max()
             msg = 'Derivative test failed at d_delay_d_%s with max relative difference %lf' % (p, np.nanmax(relative_diff).value)
             if p in ['PMELONG', 'ELONG']:
                 tol = 2e-2
             elif p in ['FB2', 'FB3']:
                 tol = 0.08
             else:
                 tol = 1e-3
             log.debug( "derivative relative diff for %s, %lf"%('d_delay_d_'+p, np.nanmax(relative_diff).value))
             assert np.nanmax(relative_diff) < tol, msg
         else:
             continue