def plot_filter(filter_type, alpha, ts, sc, overlap):
    import matplotlib.pyplot as plt
    time_taps = gfdm_filter_taps(filter_type, alpha, ts, sc)
    freq_taps = gfdm_freq_taps(time_taps)
    freq_taps_sparse = gfdm_freq_taps_sparse(freq_taps, ts, overlap)

    fig = plt.figure()
    tp = fig.add_subplot('211')
    t = np.arange(0, ts, 1. / sc)
    time_taps = np.fft.ifftshift(time_taps)
    plt.plot(t, np.abs(time_taps))
    plt.plot(t, np.abs(np.fft.ifftshift(freq_tapered_raised_cosine(t - 1. * ts / 2., alpha))))
    plt.plot(t, np.roll(np.abs(time_taps), sc))
    plt.xlim((0, ts))
    plt.xlabel('timeslot')
    plt.grid()

    fp = fig.add_subplot('212')
    f = np.arange(0, sc, 1. / ts)
    plt.plot(f, np.abs(freq_taps))
    plt.plot(f, np.abs(np.fft.fft(freq_tapered_raised_cosine(t - 1. * ts / 2., alpha))))
    plt.plot(f, np.abs(np.concatenate((freq_taps_sparse[0:len(freq_taps_sparse) / 2],
                                       np.zeros(sc * ts - len(freq_taps_sparse)),
                                       freq_taps_sparse[len(freq_taps_sparse) / 2:]))), linestyle='--')

    plt.plot(f, np.roll(np.abs(freq_taps), ts * (sc // 2)))
    plt.plot(f, np.roll(np.abs(freq_taps), ts * (sc // 2 + 1)))
    plt.xlim((0, sc))
    plt.xlabel('subcarrier')
    plt.grid()

    plt.gcf().suptitle("GFDM filter: type='{}' with M={}, K={}, L={}".format(filter_type.upper(), ts, sc, overlap), fontsize=16)
    plt.show()

def laplacian(grid, out):
    np.copyto(out, grid)
    out *= -4
    out += np.roll(grid, +1, 0)
    out += np.roll(grid, -1, 0)
    out += np.roll(grid, +1, 1)
    out += np.roll(grid, -1, 1)

 def _make_segment(self,x,y,threshold=None):
     if threshold is None:
         threshold = self._segment_threshold
     x,y=np.atleast_1d(x),np.atleast_1d(y)
     d2 = np.sqrt((np.roll(x,1)-x)**2+(np.roll(y,1)-y)**2)
     w=np.where(d2 > threshold)[0]
     #w=w[w!=0]
     xx=[]
     yy=[]
     if len(w) == 1:
         x=np.roll(x,-w[0])
         y=np.roll(y,-w[0])
         xx.append(x)
         yy.append(y)
     elif len(w) >= 2:
         xx.append(x[0:w[0]])
         yy.append(y[0:w[0]])
         for i in xrange(len(w)-1):
             xx.append(x[w[i]:w[i+1]])
             yy.append(y[w[i]:w[i+1]])
         xx.append(x[w[-1]:])
         yy.append(y[w[-1]:])
     else:
         xx.append(x)
         yy.append(y)
     return xx,yy

def make_step(net, step_size=1.5, end=default_layer, jitter=32, clip=True, objective=objective_L2, sigma=0):
    '''Basic gradient ascent step.'''

    src = net.blobs['data'] # input image is stored in Net's 'data' blob
    dst = net.blobs[end]    # the layer targeted by default_layer

    ox, oy = np.random.randint(-jitter, jitter+1, 2)
    src.data[0] = np.roll(np.roll(src.data[0], ox, -1), oy, -2) # apply jitter shift

    net.forward(end=end)     # inference of features

    objective(dst)           # set an objective

    net.backward(start=end)  # retrain

    g = src.diff[0]

    # apply normalized ascent step to the input image
    src.data[:] += step_size/np.abs(g).mean() * g

    src.data[0] = np.roll(np.roll(src.data[0], -ox, -1), -oy, -2) # unshift image

    if clip:
        bias = net.transformer.mean['data']
        src.data[:] = np.clip(src.data, -bias, 255-bias)

    if sigma:
        src.data[0] = blur(src.data[0], sigma)

def _decorate_contour_segment(data, stride=1, options={}, tomax=True, labelled=False, outline=None, aspect=1):
	default_options = {'scale': 0.2,
			'scale_units': 'dots',
			'headaxislength': 2,
			'headlength': 2,
			'headwidth': 2,
			'minshaft': 1,
			'units': 'dots',
			#'angles': 'xy',
			'edgecolor': outline,
			'linewidth': 0 if outline is None else 0.2
		}
	default_options.update(options)

	x = data[::stride,0]
	y = data[::stride,1]

	sign = 1 if tomax else -1
	dx = -sign*np.diff(y)*aspect
	dy = sign*np.diff(x)
	l = np.sqrt(dx**2+dy**2)
	dx /= l
	dy /= l

	x = 0.5*(x+np.roll(x,-1))
	y = 0.5*(y+np.roll(y,-1))

	if labelled:
		x,y,dx,dy = x[1:-2], y[1:-2], dx[1:-1], dy[1:-1]
	else:
		x,y = x[:-1], y[:-1]

	plt.quiver(x, y, dx, dy, **default_options)

    def _process_data_gated(self):
        """
        Processes the raw data from the counting device
        @return:
        """
        # remember the new count data in circular array
        self.countdata[0] = np.average(self.rawdata[0])
        # move the array to the left to make space for the new data
        self.countdata = np.roll(self.countdata, -1)
        # also move the smoothing array
        self.countdata_smoothed = np.roll(self.countdata_smoothed, -1)
        # calculate the median and save it
        self.countdata_smoothed[-int(self._smooth_window_length / 2) - 1:] = np.median(
            self.countdata[-self._smooth_window_length:])

        # save the data if necessary
        if self._saving:
            # if oversampling is necessary
            if self._counting_samples > 1:
                self._sampling_data = np.empty((self._counting_samples, 2))
                self._sampling_data[:, 0] = time.time() - self._saving_start_time
                self._sampling_data[:, 1] = self.rawdata[0]
                self._data_to_save.extend(list(self._sampling_data))
            # if we don't want to use oversampling
            else:
                # append tuple to data stream (timestamp, average counts)
                self._data_to_save.append(np.array((time.time() - self._saving_start_time,
                                                    self.countdata[-1])))
        return

def local_maxima(array2d,user_peak,index=False,count=4,floor=0,bug=False):

    from operator import itemgetter, attrgetter
    
    if user_peak == 0:
        where = ((array2d >= np.roll(array2d,1,0)) &
                 (array2d >= np.roll(array2d,-1,0)) &
                 (array2d >= np.roll(array2d,0,1)) &
                 (array2d >= np.roll(array2d,0,-1)) &
                 (array2d >= array2d.max()/5.0) &
                 (array2d > floor*np.ones(array2d.shape)) &
                 (array2d >= array2d.mean()))
    else: #some simpler filter if user indicated some modes
        where = array2d > floor

    #ignore the lesser local maxima, throw out anything with a ZERO
    if bug==True:    
        print array2d,array2d[where.nonzero()],where.nonzero()[0]
    
    peaks = zip(where.nonzero()[0],where.nonzero()[1],array2d[where.nonzero()])
    
    peaks = sorted(peaks,key=itemgetter(2),reverse=True)
   
    if len(peaks) > count and user_peak==0:
        peaks = peaks[0:count]
    
            
    keys = ['y_i','z_i','amp']
    
    peaks = [dict(zip(keys,peaks[x])) for x in range(len(peaks))]
    
    return peaks

 def GetPerimeter(self,sites=None):
     """
     Find out the associated perimeter of this cluster.
     """
     if sites is None:
         indices = self.sites.get_indices()
     else:
         indices = sites.get_indices()
     area = len(indices)
     if area <= 3:
         return 2.*(area+1.)
     else:
         cluster = numpy.zeros((self.N,self.N))
         cluster[indices[:,0],indices[:,1]] = 1.
         mask = (cluster>0)
         newcluster = numpy.copy(cluster)
         newcluster += 1.*(cluster==numpy.roll(cluster,1,0))*mask
         newcluster += 1.*(cluster==numpy.roll(cluster,1,1))*mask
         newcluster += 1.*(cluster==numpy.roll(cluster,-1,0))*mask
         newcluster += 1.*(cluster==numpy.roll(cluster,-1,1))*mask
         """
         Any site with n (1<=n<=4) neighbors will contribute (4-n) to the perimeter;
         """
         perimeter = numpy.sum(newcluster==2.)*3.+numpy.sum(newcluster==3.)*2.+numpy.sum(newcluster==4.)
         return perimeter 

	def gradient_ascent_step(net, step_size=1.5, end="inception_4c/output",
		jitter=32, clip=True, objective_fn=None, **objective_params):
		# if the objective function is None, initialize it as
		# the standard L2 objective
		if objective_fn is None:
			objective_fn = BatCountry.L2_objective

		# input image is stored in Net's 'data' blob
		src = net.blobs["data"]
		dst = net.blobs[end]

		# apply jitter shift
		ox, oy = np.random.randint(-jitter, jitter + 1, 2)
		src.data[0] = np.roll(np.roll(src.data[0], ox, -1), oy, -2)

		net.forward(end=end)
		objective_fn(dst, **objective_params)
		net.backward(start=end)
		g = src.diff[0]

		# apply normalized ascent step to the input image
		src.data[:] += step_size / np.abs(g).mean() * g

 		# unshift image
		src.data[0] = np.roll(np.roll(src.data[0], -ox, -1), -oy, -2)

		# unshift image
		if clip:
			bias = net.transformer.mean["data"]
			src.data[:] = np.clip(src.data, -bias, 255 - bias)

	def __init__(self,*arg,**kw):
		super(build_model, self).__init__(*arg, **kw)
		self['nlayers'] = 5
		self['nx'] = 500
		fault_throw = 20
		
		self['dz'] = np.array([40, 80, 40, 200, 400, ])
		self['vp'] = np.array([800., 2200., 1800., 2400., 4500., ])
		self['vs'] = self['vp']/2.
		self['rho'] = np.array([1500., 2500., 1400., 2700., 4500., ])
		self['depths'] = np.cumsum(self['dz'])
		
		self['model'] = {}
		for model in ['vp', 'vs', 'rho']:
			layer_list = []
			for index in range(self['nlayers']):
				layer = np.ones((self['nx'], self['dz'][index]), 'f')
				layer *= self[model][index]
				layer_list.append(layer)
			self['model'][model] = np.hstack(layer_list)
			self['model'][model][250:500,120:160] = self[model][1]
			self['model'][model][250:500,120+fault_throw:160+fault_throw] = self[model][2]
		
		self['model']['z'] = self['model']['vp'] * self['model']['rho']
		self['model']['R'] = (np.roll(self['model']['z'], shift=-1) - self['model']['z'])/(np.roll(self['model']['z'], shift=-1) + self['model']['z'])
		self['model']['R'][:,0] *= 0
		self['model']['R'][:,-1:] *= 0
		self['model']['R'][:,:self['dz'][0]+2] = 0

def normals_numpy(depth, rect=((0,0),(640,480)), win=7, mat=None):
    assert depth.dtype == np.float32
    from scipy.ndimage.filters import uniform_filter
    (l,t),(r,b) = rect
    v,u = np.mgrid[t:b,l:r]
    depth = depth[v,u]
    depth[depth==0] = -1e8  # 2047
    depth = calibkinect.recip_depth_openni(depth)
    depth = uniform_filter(depth, win)
    global duniform
    duniform = depth

    dx = (np.roll(depth,-1,1) - np.roll(depth,1,1))/2
    dy = (np.roll(depth,-1,0) - np.roll(depth,1,0))/2
    #dx,dy = np.array(depth),np.array(depth)
    #speedup.gradient(depth.ctypes.data, dx.ctypes.data,
    # dy.ctypes.data, depth.shape[0], depth.shape[1])

    X,Y,Z,W = -dx, -dy, 0*dy+1, -(-dx*u + -dy*v + depth).astype(np.float32)

    mat = calibkinect.projection().astype('f').transpose()
    mat = np.ascontiguousarray(mat)
    x = X*mat[0,0] + Y*mat[0,1] + Z*mat[0,2] + W*mat[0,3]
    y = X*mat[1,0] + Y*mat[1,1] + Z*mat[1,2] + W*mat[1,3]
    z = X*mat[2,0] + Y*mat[2,1] + Z*mat[2,2] + W*mat[2,3]
    w = np.sqrt(x*x + y*y + z*z)
    w[z<0] *= -1
    weights = z*0+1
    weights[depth<-1000] = 0
    weights[(z/w)<.1] = 0
    #return x/w, y/w, z/w
    return np.dstack((x/w,y/w,z/w)), weights

def test_shift():

    amp = 1
    v0 = 0 * u.m / u.s
    sigma = 8
    spectral_axis = np.arange(-50, 51) * u.m / u.s

    true_spectrum = gaussian(spectral_axis.value,
                             amp, v0.value, sigma)

    # Shift is an integer, so rolling is equivalent
    rolled_spectrum = np.roll(true_spectrum, 10)

    shift_spectrum = fourier_shift(true_spectrum, 10)

    np.testing.assert_allclose(shift_spectrum,
                               rolled_spectrum,
                               rtol=1e-4)

    # With part masked
    masked_spectrum = true_spectrum.copy()
    mask = np.abs(spectral_axis.value) <= 30
    masked_spectrum[~mask] = np.NaN

    rolled_mask = np.roll(mask, 10)
    rolled_masked_spectrum = rolled_spectrum.copy()
    rolled_masked_spectrum[~rolled_mask] = np.NaN

    shift_spectrum = fourier_shift(masked_spectrum, 10)

    np.testing.assert_allclose(shift_spectrum,
                               rolled_masked_spectrum,
                               rtol=1e-4)

def face_ibug_68_mirrored_to_face_ibug_68(pcloud):
    r"""
    Apply the IBUG 68-point semantic labels, on a pointcloud that has been
    mirrored around the vertical axis (flipped around the Y-axis). Thus, on
    the flipped image the jaw etc would be the wrong way around. This
    rectifies that and returns a new PointCloud whereby all the points
    are oriented correctly.

    The semantic labels applied are as follows:

      - jaw
      - left_eyebrow
      - right_eyebrow
      - nose
      - left_eye
      - right_eye
      - mouth

    References
    ----------
    .. [1] http://www.multipie.org/
    .. [2] http://ibug.doc.ic.ac.uk/resources/300-W/
    """
    new_pcloud, old_map = face_ibug_68_to_face_ibug_68(pcloud,
                                                       return_mapping=True)
    lms_map = np.hstack([old_map['jaw'][::-1],
                         old_map['right_eyebrow'][::-1],
                         old_map['left_eyebrow'][::-1],
                         old_map['nose'][:4],
                         old_map['nose'][4:][::-1],
                         np.roll(old_map['right_eye'][::-1], 4),
                         np.roll(old_map['left_eye'][::-1], 4),
                         np.roll(old_map['mouth'][:12][::-1], 7),
                         np.roll(old_map['mouth'][12:][::-1], 5)])
    return new_pcloud.from_vector(pcloud.points[lms_map]), old_map

    def _calculate_series(self, attr):
        mi = getattr(self, '{}_min'.format(attr))
        amp = getattr(self, '{}_max'.format(attr)) - mi

        funcname = getattr(self, '{}_func'.format(attr))
        period = getattr(self, '{}_period'.format(attr))
        offset = getattr(self, '{}_offset'.format(attr))
        speriod = getattr(self, '{}_sample'.format(attr))
        duty = getattr(self, '{}_duty'.format(attr))
        duration = getattr(self, '{}_duration'.format(attr))

        x, y, sx, sy = [], [], [], []
        if funcname != NULL_STR:
            if self.xy_pattern_enabled and getattr(self, '{}_use_transit_time'.format(attr)):
                t = self.calculate_transit_time()
            else:
                t = duration

            t = t or 1
            x = linspace(0, t, 500)

            speriod = speriod or 1
            sx = arange(0, t, speriod)
            if sx[-1] < t:
                sx = append(sx, t)

            if funcname == 'Sine':
                def func(xx):
                    return (mi + amp) + amp * sin(period * xx + offset)

                y = func(x)

            elif funcname == 'Square':
                def func(xx):
                    return 0.5 * amp * (signal.square(period * xx * 2 * pi + offset, duty=duty / 100.) + 1) + mi

                y = func(x)

                bx = asarray(diff(y), dtype=bool)
                bx = roll(bx, 1)

                sx = x[bx]
                sx = append(asarray([0]), sx)
                sx = append(sx, asarray([t]))

            elif funcname == 'Saw':
                def func(xx):
                    return 0.5 * amp * (signal.sawtooth(period * xx * 2 * pi + offset) + 1) + mi

                y = func(x)

                asign = sign(gradient(y))
                signchange = ((roll(asign, 1) - asign) != 0).astype(bool)
                signchange[0] = False
                nsx = x[signchange]
                sx = hstack((sx, nsx))

            sy = func(sx)

        return x, y, sx, sy

    def polyhedra(self, wm):
        '''Iterates through the polyhedra that make up the closest volume to a certain vertex'''
        for p, facerow in enumerate(self.connected):
            faces = facerow.indices
            pts, polys = _ptset(), _quadset()
            if len(faces) > 0:
                poly = np.roll(self.polys[faces[0]], -np.nonzero(self.polys[faces[0]] == p)[0][0])
                assert pts[wm[p]] == 0
                assert pts[self.pts[p]] == 1
                pts[wm[poly[[0, 1]]].mean(0)]
                pts[self.pts[poly[[0, 1]]].mean(0)]

                for face in faces:
                    poly = np.roll(self.polys[face], -np.nonzero(self.polys[face] == p)[0][0])
                    a = pts[wm[poly].mean(0)]
                    b = pts[self.pts[poly].mean(0)]
                    c = pts[wm[poly[[0, 2]]].mean(0)]
                    d = pts[self.pts[poly[[0, 2]]].mean(0)]
                    e = pts[wm[poly[[0, 1]]].mean(0)]
                    f = pts[self.pts[poly[[0, 1]]].mean(0)]

                    polys((0, c, a, e))
                    polys((1, f, b, d))
                    polys((1, d, c, 0))
                    polys((1, 0, e, f))
                    polys((f, e, a, b))
                    polys((d, b, a, c))

            yield pts.points, np.array(list(polys.triangles))

def get_corr_map_photon(coo1, coo2, skypos, skyrange, sec, bound, bandwidth):
  imsz = imagetools.deg2pix(skypos, skyrange, 0.0001)
  count = np.zeros(imsz)
  print(imsz)
  co_rel = np.array([[0,0]])
  len1 = coo1.shape[0]
  len2 = coo2.shape[0]
  print(len1,len2)
  wcs = imagetools.define_wcs(skypos,skyrange,width=False,height=False,verbose=0,pixsz=0.0001)

  if len2<=80 or len1<=80:
    return count, np.array([0.0, 0.0])
  if len2>len1:
    for i in range(len2):
      #print(i)
      tmp_co = np.roll(coo2, i, axis=0)[0:len1,:]-coo1
      tmp_co = tmp_co[np.absolute(tmp_co[:,0])<=bound,:]
      tmp_co = tmp_co[np.absolute(tmp_co[:,1])<=bound,:]
      co_rel = np.concatenate((co_rel, tmp_co), axis = 0)

  else:
    for i in range(len1):
      #print(i)
      tmp_co = coo2-np.roll(coo1, i, axis=0)[0:len2,:]
      tmp_co = tmp_co[np.absolute(tmp_co[:,0])<=bound,:]
      tmp_co = tmp_co[np.absolute(tmp_co[:,1])<=bound,:]
      co_rel = np.concatenate((co_rel, tmp_co), axis = 0)

  print co_rel.shape
  if co_rel.shape[0]>1000:
    centroid = ck.find_centroid(co_rel[1:], bandwidth, 11, bound)
  else:
    return count, np.array([0.0, 0.0])

  return count, centroid

	def GetReconstructionWithBoundary(self,N,orientationmap,filename):
		"""
		Decorate the orientation map with boundaries.
		"""
		indexmap = self.GetIndexmap()
		newindexmap = numpy.empty(numpy.array(indexmap.shape)*4)
		for i in range(4):
			for j in range(4):
				newindexmap[i::4,j::4] = indexmap
		neworientationmap = numpy.empty([4*N,4*N,3])
		for k in range(3):
			for i in range(4):
				for j in range(4):
					neworientationmap[i::4,j::4,k] = orientationmap[:,:,k]
		boundarymap =  1-(indexmap==numpy.roll(indexmap, 1,0))*(indexmap==numpy.roll(indexmap,-1,0))*\
					   (indexmap==numpy.roll(indexmap, 1,1))*(indexmap==numpy.roll(indexmap,-1,1))
		xs,ys = boundarymap.nonzero()
		for (i,j) in zip(xs,ys):
			temp = [indexmap[i,j]==indexmap[(i-1+N)%N,j],indexmap[i,j]==indexmap[(i+1)%N,j],\
					indexmap[i,j]==indexmap[i,(j-1+N)%N],indexmap[i,j]==indexmap[i,(j+1)%N]]
			pos = [[(4*i,4*j),(4*i+1,4*j+4)],[(4*i+3,4*j),(4*i+4,4*j+4)],\
				   [(4*i,4*j),(4*i+4,4*j+1)],[(4*i,4*j+3),(4*i+4,4*j+4)]]
			for n in range(4):
				if not temp[n]:
					newindexmap[pos[n][0][0]:pos[n][1][0],pos[n][0][1]:pos[n][1][1]] = -1
		for k in range(3):
			neworientationmap[:,:,k] *= (newindexmap!=-1)
		"""
		Use PIL to plot for larger sizes, since we want to be able to draw pixel by pixel.
		"""
		from PIL import Image
		pilImage = Image.fromarray(neworientationmap.astype('uint8'), 'RGB')
		pilImage.save(filename+'.png')

def main():
    args = parse_args()

    gdb = genome.db.GenomeDB(assembly=args.assembly)

    fwd_track = gdb.open_track(args.fwd_track)
    rev_track = gdb.open_track(args.rev_track)
    combined_track = gdb.create_track(args.combined_track)
    
    shift = find_strand_shift(gdb, fwd_track, rev_track)

    sys.stderr.write("shifting fwd/rev by +%d/-%d bp\n" % (shift, shift))

    for chrom in gdb.get_chromosomes():
        sys.stderr.write("%s\n" % chrom.name)
        carray = create_carray(combined_track, chrom, args.dtype)

        fwd_vals = fwd_track.get_nparray(chrom)
        rev_vals = rev_track.get_nparray(chrom)

        # shift fwd / rev values by the offset that gave
        # the maximum covariance2
        fwd_vals = np.roll(fwd_vals, shift)
        rev_vals = np.roll(rev_vals, -shift)
        fwd_vals[:shift] = 0
        rev_vals[-shift:] = 0

        carray[:] = fwd_vals + rev_vals
        
    fwd_track.close()
    rev_track.close()
    combined_track.close()

    def draw_svd(self,comp='V',nmodes=6,center=False,zoom = False,
                 interp='none'):
        #compute svd is not present . . later

        fig = plt.figure()
        try:
            y = self.svd[comp]
        except:
            y = self.svd['V']

        for i in  xrange(nmodes):
            canvas = fig.add_subplot(2,3,i+1)
            
            # if comp=='V':
            #     y_i = y[i,:,:]
            # else:
            #     y_i = y[:,i]
            y_i = y[i]
        
            nx,ny = y_i.shape
            if center:
                y_i = np.roll(np.roll(y_i,nx/2,axis=0),ny/2,axis=1)
            if zoom and not center:
                y_i = y_i[0:nx/8.0,0:ny/8.]
            if zoom and center:
                y_i = y_i[nx/2-nx/8:nx/2.0+nx/8,
                          ny/2-ny/16:ny/2.0+ny/16,]
    
            canvas.imshow(np.real(y_i),aspect='auto',interpolation= interp)
        
        plt.close(fig)
        return fig

def find_max_cov(fwd_vals, rev_vals):    
    max_cov = 0.0
    max_cov_shift = 0

    for i in range(MIN_SHIFT, MAX_SHIFT):
        # shift fwd values ahead by i bp
        # shift rev values back by i bp
        matrix = np.vstack((np.roll(fwd_vals, i),
                            np.roll(rev_vals, -i)))

        # compute covariance between fwd and rev
        cov = np.cov(matrix)[0,1]

        max_str = ""
        if cov > max_cov:
            # this is the highest covariance we've seen yet
            max_cov = cov
            max_cov_shift = i
            max_str = " *"

        last_max = i - max_cov_shift

        if last_max > MAX_SHIFT_MAX_DIST:
            # we seem to be well past the peak covariance
            break

        sys.stderr.write("shift: %d, cov: %g%s\n" % (i, cov, max_str))
        
    sys.stderr.write("shift: %d, cov: %g\n" % (max_cov_shift, max_cov))

    return max_cov_shift