def plotCurves(c1, c2):
    name1, t, avg1, top1, bottom1 = c1
    name2, t, avg2, top2, bottom2 = c2
    pl.plot(t, np.zeros(len(t)), 'k-')
    s1 = ma.array(avg1)
    s2 = ma.array(avg2)
    zx1 = np.logical_and(np.greater_equal(top1, 0), np.less_equal(bottom1, 0))
    zx2 = np.logical_and(np.greater_equal(top2, 0), np.less_equal(bottom2, 0))
    ix = np.logical_or(
            np.logical_and(
                np.greater_equal(top1, top2),
                np.less_equal(bottom1, top2)),
            np.logical_and(
                np.greater_equal(top1, bottom2),
                np.less_equal(bottom1, bottom2)))
    mask1 = np.logical_or(zx1, ix)
    mask2 = np.logical_or(zx2, ix)

    print mask1
    print mask2
    print zx1
    print zx2
    print ix

    pl.plot(t, s1, "k--", linewidth=1)
    pl.plot(t, s2, "k-", linewidth=1)
    s1.mask = ix
    s2.mask = ix
    pl.plot(t, s1, "k--", linewidth=3, label=name1)
    pl.plot(t, s2, "k-", linewidth=3, label=name2)
    pl.xlabel('Time (secs)')
    pl.ylabel("Pearson correlation")

   def eval(self, band, times, z=0, k=1):
      '''Evaluate, using a spline, the value of the template at specific
      times, optionally with a redshift (in the sense that the times should
      be blueshifted before interpolating.  Also returns a mask indicating
      the interpolated points (1) and the extrapolated points (0)'''
      if len(num.shape(times)) == 0:
         evt = num.array([times/(1+z)])
         scalar = 1
      else:
         evt = times/(1+z)
         scalar = 0
      if band not in self.__dict__ and band not in ['J','H','K']:
         raise AttributeError, "Sorry, band %s is not supported by dm15temp" % \
               band
      s = dm152s(self.dm15)
      if band == 'J':
         return(0.080 + evt/s*0.05104699 + 0.007064257*(evt/s)**2 - 0.000257906*(evt/s)**3,
               0.0*evt/s + 0.06, num.greater_equal(evt/s, -12)*num.less_equal(evt/s, 10)) 
      elif band == 'H':
         return(0.050 + evt/s*0.0250923 + 0.001852107*(evt/s)**2 - 0.0003557824*(evt/s)**3,
               0.0*evt/s + 0.08, num.greater_equal(evt/s, -12)*num.less_equal(evt/s, 10)) 
      elif band == 'K':
         return(0.042 + evt/s*0.02728437+ 0.003194500*(evt/s)**2 - 0.0004139377*(evt/s)**3,
               0.0*evt/s + 0.08, num.greater_equal(evt/s, -12)*num.less_equal(evt/s, 10)) 
      evd = self.tck[band].ev(evt/self.s, evt*0+self.dm15)
      eevd = self.tck['e_'+band].ev(evt/self.s, evt*0+self.dm15)
      mask = num.greater_equal(evt/self.s, -10)*num.less_equal(evt/self.s,70)

      if scalar:
         return(evd[0], eevd[0], mask[0])
      else:
         return(evd, eevd, mask)

def rel_coron_thrupt(Pmod, ref_pos):
    # Given 2-d off-axis PSF model over a set of positions,
    # comput the FWHM throughput relative to the off-axis PSF
    # at a reference (presumably peak throughput) location.
    coron_thrupt = np.empty(Pmod.shape[:-2])

    Pref = Pmod[ref_pos]
    ref_peak = np.max(Pref)
    ref_fwhm_ind = np.greater_equal(Pref, ref_peak/2)
    ref_fwhm_sum = np.sum(Pref[ref_fwhm_ind])

    if len(coron_thrupt.shape) == 2:
        for si in range(coron_thrupt.shape[0]):
            for ti in range(coron_thrupt.shape[1]):
                P = Pmod[si,ti]
                fwhm_ind = np.greater_equal(P, np.max(P)/2)
                fwhm_sum = np.sum(P[fwhm_ind])
                coron_thrupt[si,ti] = fwhm_sum/ref_fwhm_sum

    elif len(coron_thrupt.shape) == 1:
        for ti in range(coron_thrupt.shape[0]):
            P = Pmod[ti]
            fwhm_ind = np.greater_equal(P, np.max(P)/2)
            fwhm_sum = np.sum(P[fwhm_ind])
            coron_thrupt[ti] = fwhm_sum/ref_fwhm_sum

    return coron_thrupt

def ring(x, y, height, thickness, gaussian_width):
    """
    Circular ring (annulus) with Gaussian fall-off after the solid ring-shaped region.
    """
    radius = height/2.0
    half_thickness = thickness/2.0

    distance_from_origin = np.sqrt(x**2+y**2)
    distance_outside_outer_disk = distance_from_origin - radius - half_thickness
    distance_inside_inner_disk = radius - half_thickness - distance_from_origin

    ring = 1.0-np.bitwise_xor(np.greater_equal(distance_inside_inner_disk,0.0),
                              np.greater_equal(distance_outside_outer_disk,0.0))

    sigmasq = gaussian_width*gaussian_width

    if sigmasq==0.0:
        inner_falloff = x*0.0
        outer_falloff = x*0.0
    else:
        with float_error_ignore():
            inner_falloff = np.exp(np.divide(-distance_inside_inner_disk*distance_inside_inner_disk, 2.0*sigmasq))
            outer_falloff = np.exp(np.divide(-distance_outside_outer_disk*distance_outside_outer_disk, 2.0*sigmasq))

    return np.maximum(inner_falloff,np.maximum(outer_falloff,ring))

   def __call__(self, x):
      '''Interpolate at point [x].  Returns a 3-tuple: (y, mask) where [y]
      is the interpolated point, and [mask] is a boolean array with the same
      shape as [x] and is True where interpolated and False where extrapolated'''
      if not self.setup:
         self._setup()

      if len(num.shape(x)) < 1:
         scalar = True
      else:
         scalar = False

      x = num.atleast_1d(x)
      if self.realization:
         evm = num.atleast_1d(splev(x, self.realization))
         mask = num.greater_equal(x, self.realization[0][0])*\
                num.less_equal(x,self.realization[0][-1])
      else:
         evm = num.atleast_1d(splev(x, self.tck))
         mask = num.greater_equal(x, self.tck[0][0])*num.less_equal(x,self.tck[0][-1])

      if scalar:
         return evm[0],mask[0]
      else:
         return evm,mask

def curveScore(l, curve):
    lb = l["TextBB"]
    elbl = [lb[0] - 20, lb[1], lb[0], lb[3]]
    elbr = [lb[0], lb[1], lb[2] + 20, lb[3]]
    cindex = int(curve[0])
    cdata = curve[1]
    # img=cdata[elb[1]:elb[3],elb[0]:elb[2]]
    imgl = cdata[elbl[1] : elbl[3], elbl[0] : elbl[2]]
    imgr = cdata[elbr[1] : elbr[3], elbr[0] : elbr[2]]

    # print cindex,l['Text']
    # show_img(imgl)
    # show_img(imgr)
    # points from the rectangle to the left and right of the legend word that are not white or black pixels.
    lnzps = np.where(np.logical_and(np.greater_equal(imgl[:, :, 1], 0.06), np.greater_equal(imgl[:, :, 2], 0.1)))
    rnzps = np.where(np.logical_and(np.greater_equal(imgr[:, :, 1], 0.06), np.greater_equal(imgr[:, :, 2], 0.1)))

    if len(lnzps[0]) == 0 and len(rnzps[0]) == 0:  # this means for this legend, we did not find a single pixel from the
        # curve that is to the left or right of it.
        # print l['Text'],"has no points to the left or right for",cindex
        return (None, None, None)
    elif len(lnzps[0]) != 0 and len(rnzps[0]) == 0:
        # print l['Text'],"has curve",cindex,"to the left of it, distance:",100-np.sort(lnzps[0])[-1]
        return (cindex, "l", 20 - np.sort(lnzps[0])[-1])
    elif len(lnzps[0]) == 0 and len(rnzps[0]) != 0:
        # print l['Text'],"has curve",cindex,"to the right of it, distance:",np.sort(rnzps[0])[0]
        return (cindex, "r", np.sort(rnzps[0])[0])
    else:  # this means, some points from this curve belongs to both left and the right of the legend. That is improbable.
        print "Something wrong, a single curve has pixels on both sides of the image "
        return (None, None)

def radial_contrast_flr(image, xc, yc, seps, zw, coron_thrupt, klip_thrupt=None):
    rad_flr_ctc = np.empty((len(seps)))
    assert(len(seps) == len(coron_thrupt))
    if klip_thrupt is not None:
        assert(len(seps) == len(klip_thrupt))
        rad_flr_ctc_ktc = np.empty((len(seps)))
    else:
        rad_flr_ctc_ktc = None

    imh = image.shape[0]
    imw = image.shape[1]

    xs = np.arange(imw) - xc
    ys = np.arange(imh) - yc
    XXs, YYs = np.meshgrid(xs, ys)
    RRs = np.sqrt(XXs**2 + YYs**2)

    for si, sep in enumerate(seps):
        r_in = np.max([seps[0], sep-zw/2.])
        r_out = np.min([seps[-1], sep+zw/2.])
        meas_ann_mask = np.logical_and(np.greater_equal(RRs, r_in),
                                          np.less_equal(RRs, r_out))
        meas_ann_ind = np.nonzero(np.logical_and(np.greater_equal(RRs, r_in).ravel(),
                                                    np.less_equal(RRs, r_out).ravel()))[0]
        meas_ann = np.ravel(image)[meas_ann_ind]
        rad_flr_ctc[si] = np.nanstd(meas_ann)/coron_thrupt[si]
        if rad_flr_ctc_ktc is not None:
            rad_flr_ctc_ktc[si] = np.nanstd(meas_ann)/coron_thrupt[si]/klip_thrupt[si]

    #pdb.set_trace()
    return rad_flr_ctc, rad_flr_ctc_ktc

def filter(mask, cube, header, clipMethod, threshold, rmsMode, verbose):
    if clipMethod == 'relative':
        # determine the clip level
		# Measure noise in original cube
		# rms = GetRMS(cube,rmsmode=rmsMode,zoomx=1,zoomy=1,zoomz=100000,verb=verbose,nrbins=100000)
        rms = GetRMS(cube, rmsMode=rmsMode, zoomx=1, zoomy=1, zoomz=1, verbose=verbose)
        print 'Estimated rms = ', rms
        clip = threshold * rms
    if clipMethod == 'absolute':
        clip = threshold
    print 'using clip threshold: ', clip
	#return ((cube >= clip)+(cube <= -1*clip))

   
    # check whether there are NaNs
    nan_mask = np.isnan(cube)
    found_nan=nan_mask.sum()
    if found_nan:
        cube=np.nan_to_num(cube)
        np.logical_or(mask, (np.greater_equal(cube, clip) + np.less_equal(cube, -clip)), mask)
        cube[nan_mask]=np.nan
    else:
        np.logical_or(mask, (np.greater_equal(cube, clip) + np.less_equal(cube, -clip)), mask)

	
    return 

def getMaxPoints(arr):
    # [TODO] Work out for RGB rather than array, and maybe we don't need the filter, but hopefully speeds it up.
    # Reference http://scipy-cookbook.readthedocs.io/items/FiltFilt.html
    arra = filtfilt(b,a,arr)
    maxp = maxpoints(arra, order=(len(arra)/20), mode='wrap')
    minp = minpoints(arra, order=(len(arra)/20), mode='wrap')

    points = []

    for i in range(3):
        mas = np.equal(np.greater_equal(maxp,(i*(len(arra)/3))), np.less_equal(maxp,((i+1)*len(arra)/3)))
        k = np.compress(mas[0], maxp)
        if len(k)==0:
            continue
        points.append(sum(k)/len(k))

    if len(points) == 1:
        return points, []

    points = np.compress(np.greater_equal(arra[points],(max(arra)-min(arra))*0.40 + min(arra)),points)
    rifts = []
    for i in range(len(points)-1):
        mas = np.equal(np.greater_equal(minp, points[i]),np.less_equal(minp,points[i+1]))
        k = np.compress(mas[0], minp)
        rifts.append(k[arra[k].argmin()])

    return points, rifts

def _calc_uncorr_gene_score(gene, input_gene, input_snp, pruned_snps, hotspots):
    # find local snps given a gene
    cond_snps_near_gene = logical_and(np.equal(input_snp[:, 0], input_gene[gene, 0]),
                                      np.greater_equal(input_snp[:, 1], input_gene[gene, 1]),
                                      np.less_equal(input_snp[:, 1], input_gene[gene, 2]))
    # if no snps found
    if not np.any(cond_snps_near_gene):
        return (np.nan, 0, 1, 0, 0)

    n_snps_zscore_finite = np.sum(np.isfinite(input_snp[cond_snps_near_gene][:, 3]))
    # if no snps with finite zcore
    if n_snps_zscore_finite == 0:
        return (np.nan, 0, 1, 0, 0)

    n_snps_per_gene = n_snps_zscore_finite

    # use p-value to find most significant SNP
    idx_min_pval = np.nanargmin(input_snp[cond_snps_near_gene][:, 3])

    uncorr_score = input_snp[cond_snps_near_gene][idx_min_pval, 2]

    # count number of independent SNPs per gene
    n_indep_snps_per_gene = np.sum(logical_and(np.equal(pruned_snps[:, 0], input_gene[gene, 0]),
                                               np.greater_equal(pruned_snps[:, 1], input_gene[gene, 1]),
                                               np.less_equal(pruned_snps[:, 1], input_gene[gene, 2])))

    # count number of hotspots per gene
    n_hotspots_per_gene = np.sum(np.logical_and(np.equal(hotspots[:, 0], input_gene[gene, 0]),
                                                np.greater(np.fmin(hotspots[:, 2], input_gene[gene, 2])
                                                           - np.fmax(hotspots[:, 1], input_gene[gene, 1]), 0)))
    return (uncorr_score, n_snps_per_gene, 0, n_indep_snps_per_gene, n_hotspots_per_gene)

def minima_in_range(r, g_r, r_min, r_max):
    """Find the minima in a range of r, g_r values"""
    idx = np.where(np.logical_and(np.greater_equal(r, r_min), np.greater_equal(r_max, r)))
    g_r_slice = g_r[idx]
    g_r_min = g_r_slice[g_r_slice.argmin()]
    idx_min, _ = find_nearest(g_r, g_r_min)
    return r[idx_min], g_r[idx_min]

def parallel_point_test(center,dim,x,y,z):
    '''
    Overview:
        Determines whether a given point is in a parallelapiped given the point
    being tested and the relevant parameters.


    Parameters:

    center:(float,[3]|angstroms) = The coordinates of the center of the
    parallelapiped. This parameter is in the form (x center,y center, z center)

    dim:(float,[3]|angstroms) = The x, y and z dimensions of the parallelapiped
    object.

    x,y,z:(float|angstroms) = coordinates for the point being tested.


    Note:
    -The API is left intentionally independent of the class structures used in
    sample_prep.py to allow for code resuabilitiy.

    '''

    low_lim = (array(center) - (array(dim)/2.0))
    high_lim = (array(center) +(array(dim)/2.0))

    height_lim = greater_equal (z,low_lim[2])*less_equal (z,high_lim[2])
    length_lim = greater_equal (y,low_lim[1])*less_equal (y,high_lim[1])
    width_lim = greater_equal (x,low_lim[0])*less_equal (x,high_lim[0])

    test_results = height_lim * length_lim * width_lim

    return test_results

def arc_by_radian(x, y, height, radian_range, thickness, gaussian_width):
    """
    Radial arc with Gaussian fall-off after the solid ring-shaped
    region with the given thickness, with shape specified by the
    (start,end) radian_range.
    """

    # Create a circular ring (copied from the ring function)
    radius = height/2.0
    half_thickness = thickness/2.0

    distance_from_origin = np.sqrt(x**2+y**2)
    distance_outside_outer_disk = distance_from_origin - radius - half_thickness
    distance_inside_inner_disk = radius - half_thickness - distance_from_origin

    ring = 1.0-np.bitwise_xor(np.greater_equal(distance_inside_inner_disk,0.0),
                              np.greater_equal(distance_outside_outer_disk,0.0))

    sigmasq = gaussian_width*gaussian_width

    if sigmasq==0.0:
        inner_falloff = x*0.0
        outer_falloff = x*0.0
    else:
        with float_error_ignore():
            inner_falloff = np.exp(np.divide(-distance_inside_inner_disk*distance_inside_inner_disk, 2.0*sigmasq))
            outer_falloff = np.exp(np.divide(-distance_outside_outer_disk*distance_outside_outer_disk, 2.0*sigmasq))

    output_ring = np.maximum(inner_falloff,np.maximum(outer_falloff,ring))

    # Calculate radians (in 4 phases) and cut according to the set range)

    # RZHACKALERT:
    # Function float_error_ignore() cannot catch the exception when
    # both np.dividend and divisor are 0.0, and when only divisor is 0.0
    # it returns 'Inf' rather than 0.0. In x, y and
    # distance_from_origin, only one point in distance_from_origin can
    # be 0.0 (circle center) and in this point x and y must be 0.0 as
    # well. So here is a hack to avoid the 'invalid value encountered
    # in divide' error by turning 0.0 to 1e-5 in distance_from_origin.
    distance_from_origin += np.where(distance_from_origin == 0.0, 1e-5, 0)

    with float_error_ignore():
        sines = np.divide(y, distance_from_origin)
        cosines = np.divide(x, distance_from_origin)
        arcsines = np.arcsin(sines)

    phase_1 = np.where(np.logical_and(sines >= 0, cosines >= 0), 2*pi-arcsines, 0)
    phase_2 = np.where(np.logical_and(sines >= 0, cosines <  0), pi+arcsines,   0)
    phase_3 = np.where(np.logical_and(sines <  0, cosines <  0), pi+arcsines,   0)
    phase_4 = np.where(np.logical_and(sines <  0, cosines >= 0), -arcsines,     0)
    arcsines = phase_1 + phase_2 + phase_3 + phase_4

    if radian_range[0] <= radian_range[1]:
        return np.where(np.logical_and(arcsines >= radian_range[0], arcsines <= radian_range[1]),
                        output_ring, 0.0)
    else:
        return np.where(np.logical_or(arcsines >= radian_range[0], arcsines <= radian_range[1]),
                        output_ring, 0.0)

def cone_point_test(center,dim,stub,x,y,z):
    '''
    Overview:
        Determines whether a given point is in an cone given the point being
    tested and the relevant parameters..


    Parameters:

    center:float,[3]|angstroms) = The x, y, and z component of the central
    point of the ellipsoid. In the case that the center is set to
    [None,None,None] the shape will be put in the bottom corner of the unit cell
    (the bounding box will start at (0,0,0).

    dim:(float,[3]|angstroms) = The x component, y component and thickness
    of the cone respectively. x is the radius of the cone base in the x
    direction and b is the radius of the cone base in the y direction.

    stub:(float|angstroms) = provides a hard cut-off for the thickness of the
    cone. this allows for the creation of a truncated cone object who side slope
    can be altered by using different z component values while keeping the stub
    parameter fixed.

    x,y,z:(float|angstroms) = coordinates for the point being tested.


    Notes:
    -To solve this equation more efficiently, the program takes in an array of
    x,y and z so that x[size(x),1,1], y[1,size(y),1], z[1,1,size(z)]. This
    module then solves each part of the test individually and takes the product.
    Only the points where all of the inquires are True will be left as true in
    the test_results array

    -The API is left intentionally independent of the class structures used in
    sample_prep.py to allow for code resuabilitiy.

    '''

    a_angle = arctan(dim[2]/dim[0])
    b_angle = arctan(dim[2]/dim[1])

    low_height_lim = greater_equal (z,(center[2] - dim[2]/2))

    if stub == None:
        up_height_lim =  less_equal (z,(center[2] + dim[2]/2))
    else:
        up_height_lim =  less_equal (z,(center[2] + stub/2))

    xy_test = ((((x-center[0])**2)/((((center[2] +
           dim[2]/2)-z)/tan(a_angle))**2))+(((y-center[1])**2)/((((center[2] +
           dim[2]/2)-z)/tan(b_angle))**2)))

    in_plane_low_lim = less_equal (0.0,xy_test)
    in_plane_high_lim = greater_equal (1.0,xy_test)

    test_results = (low_height_lim * up_height_lim * in_plane_low_lim *
                    in_plane_high_lim)

    return test_results

 def __ge__(a, b):
     try:
         return np.greater_equal(a.v, b.v)
     except AttributeError:
         if isinstance(a, Measurement):
             return np.greater_equal(a.v, b)
         else:
             return np.greater_equal(a, b.v)

def gradient_to_spherical(gx,gy):
	"""
	This function convert gradient coordinates of the 
	reflector into spherical coordinates of reflected rays
	on the unit sphere S2.
	
	Parameters
	----------
	gx : 1D array
		Gradients coordinate along x axis
	gy : 1D array
		Gradients coordinate along y axis
	
	Returns
	-------
	theta : 1D array
		Inclination angles (with respect to the
		positiv z axis). 0 <= theta <= pi
	phi : 1D array
		Azimuthal angles (projection of a direction
		in z=0 plane with respect to the x axis).
		0 <= phi <= 2pi
		
	See Also
	--------
	Inverse Methods for Illumination Optics, Corien Prins
	"""
	try:
		if len(gx.shape) > 1 or len(gy.shape) > 1:
			raise NotProperShapeError("gx and gy must be 1D arrays.")
		
		if gx.shape != gy.shape:
			raise NotProperShapeError("gx and gy must have the same length.")
			
		# theta computation
		num = gx*gx + gy*gy - 1
		denom = gx*gx + gy*gy + 1
		theta = np.arccos(num/denom)

		# phi computation
		zero = np.zeros(gx.shape)
		phi = np.zeros(gx.shape)
		J = np.logical_and(np.greater_equal(gx,zero),np.greater_equal(gy,zero))
		phi[J] = np.arctan(gy[J]/gx[J])
			
		J = np.less(gx, zero)
		phi[J] = np.arctan(gy[J]/gx[J]) + np.pi
			
		J = np.logical_and(np.greater_equal(gx, zero), np.less(gy, zero))
		phi[J] = np.arctan(gy[J]/gx[J]) + 2*np.pi
			
		return theta, phi
		
	except FloatingPointError:
		print("****gradient_to_spherical error: division by zero.")
		
	except NotProperShapeError, arg:
		print("****gradient_to_spherical error: ", arg.msg)

    def _numpy(self, data, weights, shape):
        q = self.quantity(data)
        self._checkNPQuantity(q, shape)
        self._checkNPWeights(weights, shape)
        weights = self._makeNPWeights(weights, shape)
        newentries = weights.sum()

        import numpy

        selection = numpy.isnan(q)
        numpy.bitwise_not(selection, selection)
        subweights = weights.copy()
        subweights[selection] = 0.0
        self.nanflow._numpy(data, subweights, shape)

        # avoid nan warning in calculations by flinging the nans elsewhere
        numpy.bitwise_not(selection, selection)
        q = numpy.array(q, dtype=numpy.float64)
        q[selection] = self.high
        weights = weights.copy()
        weights[selection] = 0.0

        numpy.greater_equal(q, self.low, selection)
        subweights[:] = weights
        subweights[selection] = 0.0
        self.underflow._numpy(data, subweights, shape)

        numpy.less(q, self.high, selection)
        subweights[:] = weights
        subweights[selection] = 0.0
        self.overflow._numpy(data, subweights, shape)

        if all(isinstance(value, Count) and value.transform is identity for value in self.values) and numpy.all(numpy.isfinite(q)) and numpy.all(numpy.isfinite(weights)):
            # Numpy defines histograms as including the upper edge of the last bin only, so drop that
            weights[q == self.high] == 0.0

            h, _ = numpy.histogram(q, self.num, (self.low, self.high), weights=weights)

            for hi, value in zip(h, self.values):
                value.fill(None, float(hi))

        else:
            q = numpy.array(q, dtype=numpy.float64)
            numpy.subtract(q, self.low, q)
            numpy.multiply(q, self.num, q)
            numpy.divide(q, self.high - self.low, q)
            numpy.floor(q, q)
            q = numpy.array(q, dtype=int)

            for index, value in enumerate(self.values):
                numpy.not_equal(q, index, selection)
                subweights[:] = weights
                subweights[selection] = 0.0
                value._numpy(data, subweights, shape)

        # no possibility of exception from here on out (for rollback)
        self.entries += float(newentries)

def ts_increments(ts, monotony = 'increasing', max_value = None, reset_value = 0.):

    '''Return a timeserie with the increments registered in the 
        input timeserie

    .. arguments:
    - (list) ts: pandas DataFrame containing a timeserie
    - (string) monotony: increasing / decreasing / non_monotonous
    - (float) max_value: value from which the meter is reseted
    - (float) reset_value: value to which the meter is reseted

    .. returns:
    - on success: timeseries of increments. The output timeseries contains 
        one value less than the original one. The diference between 
        two values, is assigned to the epoch of the second one.'''

    new_ts = ts_to_float(ts)

    if 'error' in new_ts:
        return new_ts

    if len(new_ts) <= 1:
        return {'error': 'timeserie must have length greater than 1 to compute increments'}

    if max_value != None:
        try:
            max_value = float(max_value)
        except:
            return {'error': 'max_value is not a number'}

    try:
        reset_value = float(reset_value)
    except:
        return {'error': 'reset_value is not a number'}

    if monotony == 'increasing':
        if not np.greater_equal(new_ts['value'], reset_value).all():
            return {'error': 'value lower than reset_value'}
        elif max_value and not np.less_equal(new_ts['value'], max_value).all():
            return {'error': 'value greater than max_value'}
    elif monotony == 'decreasing':
        if not np.less_equal(new_ts['value'].values, reset_value).all():
            return {'error': 'value greater than reset value'}
        elif max_value and not np.greater_equal(new_ts['value'], max_value).all():
            return {'error': 'value lower than max_value'}

    new_ts['old_value'] = new_ts['value'].shift()

    new_ts = new_ts.drop(new_ts.index[0])

    new_ts['increments'] = new_ts.apply(single_inc, axis = 1, monotony = monotony, \
        max_value = max_value, reset_value = reset_value)

    output_ts = pd.DataFrame()
    output_ts['value'] = new_ts['increments']

    return output_ts

def threshScore(reg_op, ref):
    if (reg_op.shape != ref.shape):
        raise ValueError("Scan shapes must be identical to compare.")
    threshold = 10
    binary_reg_op = np.greater_equal(reg_op,threshold)
    binary_ref = np.greater_equal(ref,threshold)
    boundary_match = np.equal(binary_reg_op,binary_ref)
    score = float(np.sum(boundary_match)) / boundary_match.size
    return score

def size_stats(l,x, msg=""):
    '''l and x are lists'''
    print msg
    # sort sizes
    l.sort()

    ## get X values for NX stats and sort
    x.sort()

    ## Get N reads
    N = len(l)
    print "N = %d" % (N)

    ## Get sum of data
    A = sum(l)
    print "Total length = %d" % (A)

    ## Get max length:
    MAX = max(l)
    print "Max length = %d" % (MAX)

    ## Get min length:
    MIN = min(l)
    print "Min length = %d" % (MIN)

    ## Get mean length
    MEAN = np.mean(l)
    print "Mean length = %d" % (MEAN)

    ## Get median contig size
    MEDIAN = np.median(l)
    print "Median length = %d" % (MEDIAN)

    ## Get NX values
    nxvalues = NX(l,x,G=A)
    for e in x:
        print "N%s length = %d" % (str(e), nxvalues[e])

    ## expected read length
    E = e_size(l,G=A)
    print "Expected size = %d" % (E)

    ##number reads >= X 
    print "N reads >= 10kb:"
    print sum(np.greater_equal(l,50e3))
    print "N reads >= 25kb:"
    print sum(np.greater_equal(l,50e3))
    print "N reads >= 50kb:"
    print sum(np.greater_equal(l,50e3))
    print "N reads >= 75kb:"
    print sum(np.greater_equal(l,75e3))
    print "N reads >= 100kb:"
    print sum(np.greater_equal(l,100e3))
    ##TODO: also want size data from reads >=X
    ## ALSO in diff fxn if both length and Q available - do longest with Q>x, etc
    print