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 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 __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 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 phantom(x):
        result = True

        for xi, xmin, xmax in zip(x, min_pt, max_pt):
            result = (result &
                      np.less_equal(xmin, xi) & np.less_equal(xi, xmax))
        return result

 def test_rand(self):
     # Simple distributional checks for sparse.rand.
     for random_state in None, 4321, np.random.RandomState():
         x = sprand(10, 20, density=0.5, dtype=np.float64,
                    random_state=random_state)
         assert_(np.all(np.less_equal(0, x.data)))
         assert_(np.all(np.less_equal(x.data, 1)))

    def phan(x):
        result = True

        for xp, minv, maxv in zip(x, begin, end):
            result = (result &
                      np.less_equal(minv, xp) & np.less_equal(xp, maxv))
        return result

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 FDR(p,q):
	'''
	input:
	p   - vector of p-values (numpy array)
	q   - False Discovery Rate level
	output:
	pID - p-value threshold based on independence or positive dependence
	pN  - Nonparametric p-value threshold
	*if non of the p-values pass FDR, an empty list is returned
	'''
	isort = np.argsort(p)
	p = p[isort]
	V = len(p)
	I = np.array(range(1,V+1))
	cVID = 1
	cVN = sum(np.divide(float(1),I))
	#threshold based on independence or positive dependence
	ID_thresh_vec = (I*q)/V/cVID
	ID_pass_vec = np.less_equal(p,ID_thresh_vec)
	if any(ID_pass_vec):
		pID = max(p[ID_pass_vec])
	else:
		pID = []
	# Nonparametric threshold
	N_thresh_vec = (I*q)/V/cVN
	N_pass_vec = np.less_equal(p,N_thresh_vec)
	if any(N_pass_vec):
		pN = max(p[N_pass_vec])
	else:
		pN = []
	return pID, pN

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 find_pi(m):
	np.random.seed()
	x1 = np.random.random(N)
	y1 = np.random.random(N)
	np.square(x1,x1)
	np.square(y1,y1)
	np.add(x1,y1,x1)
	np.less_equal(x1, 1.0, x1)
	return np.add.reduce(x1)*4.0/N

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 _subset(self, z):
     """
     Hampel's function is defined piecewise over the range of z
     """
     z = np.fabs(np.asarray(z))
     t1 = np.less_equal(z, self.a)
     t2 = np.less_equal(z, self.b) * np.greater(z, self.a)
     t3 = np.less_equal(z, self.c) * np.greater(z, self.b)
     return t1, t2, t3

def next_step(state):
    neighbors = game_tmp
    ndimage.convolve(state, kernel, output=neighbors)
    np.greater_equal(neighbors, 2, out=game_bool1)
    np.less_equal(neighbors, 3, out=game_bool2)
    np.multiply(game_bool1, game_bool2, out=game_bool1)
    np.multiply(state, game_bool1, out=state)
    np.equal(neighbors, 3, out=game_bool1)
    np.add(state, game_bool1, out=state)
    np.clip(state + game_bool1, 0, 1, out=state)

def test_universal():
    np.random.seed(98052)
    builtin_sgd = lambda params: sgd(params, lr=learning_rate_schedule(0.125, UnitType.minibatch))
    builtin_last_avg_error, builtin_avg_error, _ = ffnet(builtin_sgd)
    np.random.seed(98052)
    my_sgd = lambda ps, gs: C.combine([C.assign(p, p - 0.125/25 * g) for p, g in zip(ps, gs)])
    universal_sgd = lambda params: universal(my_sgd, params)
    my_last_avg_error, my_avg_error, _ = ffnet(universal_sgd)
    assert np.all(np.less_equal(my_last_avg_error, builtin_last_avg_error))
    assert np.all(np.less_equal(my_avg_error, builtin_avg_error))

  def test_predict_gives_correct_shapes_in_inference_mode_first_stage_only(
      self):
    test_graph = tf.Graph()
    with test_graph.as_default():
      model = self._build_model(
          is_training=False, first_stage_only=True, second_stage_batch_size=2)
      batch_size = 2
      height = 10
      width = 12
      input_image_shape = (batch_size, height, width, 3)

      preprocessed_inputs = tf.placeholder(dtype=tf.float32,
                                           shape=(batch_size, None, None, 3))
      prediction_dict = model.predict(preprocessed_inputs)

      # In inference mode, anchors are clipped to the image window, but not
      # pruned.  Since MockFasterRCNN.extract_proposal_features returns a
      # tensor with the same shape as its input, the expected number of anchors
      # is height * width * the number of anchors per location (i.e. 3x3).
      expected_num_anchors = height * width * 3 * 3
      expected_output_keys = set([
          'rpn_box_predictor_features', 'rpn_features_to_crop', 'image_shape',
          'rpn_box_encodings', 'rpn_objectness_predictions_with_background',
          'anchors'])
      expected_output_shapes = {
          'rpn_box_predictor_features': (batch_size, height, width, 512),
          'rpn_features_to_crop': (batch_size, height, width, 3),
          'rpn_box_encodings': (batch_size, expected_num_anchors, 4),
          'rpn_objectness_predictions_with_background':
          (batch_size, expected_num_anchors, 2),
          'anchors': (expected_num_anchors, 4)
      }

      init_op = tf.global_variables_initializer()
      with self.test_session() as sess:
        sess.run(init_op)
        prediction_out = sess.run(prediction_dict,
                                  feed_dict={
                                      preprocessed_inputs:
                                      np.zeros(input_image_shape)
                                  })

        self.assertEqual(set(prediction_out.keys()), expected_output_keys)

        self.assertAllEqual(prediction_out['image_shape'], input_image_shape)
        for output_key, expected_shape in expected_output_shapes.iteritems():
          self.assertAllEqual(prediction_out[output_key].shape, expected_shape)

        # Check that anchors are clipped to window.
        anchors = prediction_out['anchors']
        self.assertTrue(np.all(np.greater_equal(anchors, 0)))
        self.assertTrue(np.all(np.less_equal(anchors[:, 0], height)))
        self.assertTrue(np.all(np.less_equal(anchors[:, 1], width)))
        self.assertTrue(np.all(np.less_equal(anchors[:, 2], height)))
        self.assertTrue(np.all(np.less_equal(anchors[:, 3], width)))

def wigley(time):
  y = numpy.zeros(time.shape).astype(type(time))
  numpy.putmask(y,numpy.less(time,1953.),         \
                ((time-1860.)/(1953.-1860.)) * 35.0)
  numpy.putmask(y,numpy.logical_and(numpy.greater_equal(time,1953.),  \
                                    numpy.less_equal(time,1973.)),    \
                ((time-1953.)/(1973.-1953.)) * (68. - 35.) + 35.)
  numpy.putmask(y,numpy.logical_and(numpy.greater(time,1973.),        \
                                    numpy.less_equal(time,1990.)),    \
                ((time-1973.)/(1990.-1973.)) * (75. - 68.) + 68.)
  return y