def calc_link_dis(base_side, side1, side2):
    # print base_side.shape, side1.shape, side2.shape
    ans = np.zeros_like(base_side, dtype=np.float)
    mask = np.ones_like(base_side, dtype=np.bool)

    #point on the link
    mask_on_line = np.logical_and(base_side == side1+side2, mask)
    mask = np.logical_xor(mask, mask_on_line)
    ans[mask_on_line] = 0

    #the adjaceny points on the link is overlapped
    mask_point = np.logical_and(base_side < 1e-10, mask)
    mask = np.logical_xor(mask, mask_point)
    ans[mask_point] = side1[mask_point]

    side1_sqr = side1 * side1
    side2_sqr = side2 * side2
    base_side_sqr = base_side * base_side

    #obtuse case 1
    mask_obtuse1 = np.logical_and(side1_sqr > base_side_sqr + side2_sqr, mask)
    mask = np.logical_xor(mask, mask_obtuse1)
    ans[mask_obtuse1] = side2[mask_obtuse1]

    #obtuse case 2
    mask_obtuse2 = np.logical_and(side2_sqr > base_side_sqr + side1_sqr, mask)
    mask = np.logical_xor(mask, mask_obtuse2)
    ans[mask_obtuse2] = side1[mask_obtuse2]

    #compute height by Heron's formula
    half_p = (base_side[mask] + side1[mask] + side2[mask]) * 0.5 # half perimeter
    area = np.sqrt(half_p * (half_p - side1[mask]) * (half_p - side2[mask]) * (half_p - base_side[mask]))
    ans[mask] = 2 * area / base_side[mask]
    return ans

def generate_sky_model_alms(fits_file,lmax=10):
    # http://healpy.readthedocs.org/en/latest/generated/healpy.sphtfunc.map2alm.html#healpy.sphtfunc.map2alm
    healmap = a.map.Map(fromfits=fits_file)
    as_pos = hp.sphtfunc.map2alm(healmap.map.map, lmax=lmax, pol=False)
    alms_pos = n.zeros([as_pos.shape[0],3],dtype='complex')
    #print alms_pos.shape
    kk=0
    for ll in range(lmax+1):
        for mm in range(0,ll+1):
            alms_pos[kk] = n.array([ll,mm,as_pos[kk]])
            kk+=1
    #print alms_pos
    alms = n.zeros([(lmax+1)**2,3],dtype='complex')
    kk=0
    for ll in range(lmax+1):
        for mm in range(-ll,ll+1):
            if mm<0:
                alm = alms_pos[n.where(n.logical_and(alms_pos[:,0]==ll, alms_pos[:,1]==-mm)),2]
                #print 'less',ll,mm,alm
                alms[kk] = n.array([ll,mm,n.conj(alm[0,0])])
            else:
                alm = alms_pos[n.where(n.logical_and(alms_pos[:,0]==ll, alms_pos[:,1]==mm)),2]
                #print 'greater ',ll,mm,alm 
                alms[kk] = n.array([ll,mm,alm[0,0]])
            kk+=1
    return alms 

def events_in_polygon(locations, source_zone, flag_vector = None, 
    upper_depth=None, lower_depth=None):
    '''Function to identify valid events inside a polygon
    :param locations: xyz cartesian represent hypocentres
    :type locations: numpy.ndarray
    :param source_zone: source zone polygon in xyz format
    :type source_zone: numpy.ndarray
    
    '''
    neq = np.shape(locations)[0]
    if instance(flag_vector, np.ndarray):
        '''A flag vector is input'''
        if len(flag_vector) != neq:
            raise ValueError(
                'Flag vector length is not equal to number of events')
    else:
        '''A flag vector is needed'''
        flag_vector = np.zeros(neq, dtype = int)

    valid_id = flag_vector == 0
    if upper_depth:
        valid_id = np.logical_and(valid_id, locations[:, -1] <= upper_depth)
    
    if lower_depth:
        valid_id = np.logical_and(valid_id, locations[:, -1] >= lower_depth)
    
    valid_id[valid_id] = points_in_poly(locations[valid_id, :-1], 
                                source_zone)
    
    #if not(np.all(valid_id)):
    flag_vector[np.logical_not(valid_id)] = 1
    return flag_vector

def doCombinations(stats, varname, data, amoeboids, mesenchymals, successful, function=np.mean, axis=None):
    """Apply a given function to all elements of a given array, and do so for all desired combinations.
    The combinations are for example "all agents that are amoeboid and were successful", which would be denoted by `_a_s`.
    Note that some of these results may make little sense because too few agents fit the criteria.
    """
    a_s = np.logical_and(amoeboids, successful)
    m_s = np.logical_and(mesenchymals, successful)
    a_us = np.logical_and(amoeboids, ~successful)
    m_us = np.logical_and(mesenchymals, ~successful)
    combinations = {
                    "" : np.ones_like(data, dtype=np.bool_),
                    "_a" : amoeboids,
                    "_m" : mesenchymals,
                    "_a_s" : a_s,
                    "_m_s" : m_s,
                    "_a_us" : a_us,
                    "_m_us" : m_us,
                    }
    for suffix, selector in combinations.iteritems():
        if axis is None:
            myselector = np.s_[selector]
        elif axis==1:
            myselector = np.s_[:,selector]
        value = function( data[myselector] ) if selector.any() else None
        if value is not None and type(value)!=float:
            value = float(value)
        stats[varname + suffix] = value 
    return stats

 def evaluate(x, y, amplitude, x_0, y_0):
     """Two dimensional delta model function"""
     dx = x - x_0
     dy = y - y_0
     x_mask = np.logical_and(dx > -0.5, dx <= 0.5)
     y_mask = np.logical_and(dy > -0.5, dy <= 0.5)
     return np.select([np.logical_and(x_mask, y_mask)], [amplitude])

  def check_cross_amplification(self, max_dist):
    pop_set = set()
    for pos in self.pairs_pos:
      p1,p2 = pos

      c1,a1,s1 = zip(*self.alignments_dict[p1])
      c2,a2,s2 = zip(*self.alignments_dict[p2])
     
      chrs_x, chrs_y = na.meshgrid(c1, c2)
      chrs_eq = chrs_x==chrs_y
    
      orient_x,orient_y = na.meshgrid(s1, s2)
      orient_fr = orient_x==na.logical_not(orient_y)
    
      pos_x,pos_y = na.meshgrid(a1, a2)
    
      # Filter by minimum distance.
      # And if orientation is forward/reverse
      # Amps is [Alignments i X Alignments j ]
      amps = na.logical_and(na.logical_and(na.abs(pos_x-pos_y)<max_dist, 
                      na.logical_and(orient_fr, pos_x<pos_y)), # Forward reverse if one position is less than the other and is forward, while the other is reverse              
              chrs_eq)
      number_of_amps = na.sum(amps)
      #assert number_of_amps!=0
      if number_of_amps>1:
        pop_set.add(pos)
        
    self.pairs_pos.difference_update(pop_set)

    def set_boundary_pixels(self, value=0.0, n_pixels=1):
        r"""
        Returns a copy of this :map:`MaskedImage` for which n pixels along
        the its mask boundary have been set to a particular value. This is
        useful in situations where there is absent data in the image which
        can cause, for example, erroneous computations of gradient or features.

        Parameters
        ----------
        value : float or (n_channels, 1) ndarray
        n_pixels : int, optional
            The number of pixels along the mask boundary that will be set to 0.

        Returns
        -------
         : :map:`MaskedImage`
            The copy of the image for which the n pixels along its mask
            boundary have been set to a particular value.
        """
        global binary_erosion
        if binary_erosion is None:
            from scipy.ndimage import binary_erosion  # expensive
        # Erode the edge of the mask in by one pixel
        eroded_mask = binary_erosion(self.mask.mask, iterations=n_pixels)

        # replace the eroded mask with the diff between the two
        # masks. This is only true in the region we want to nullify.
        np.logical_and(~eroded_mask, self.mask.mask, out=eroded_mask)
        # set all the boundary pixels to a particular value
        self.pixels[..., eroded_mask] = value

def calc_India_Burma_througt(Field,I_Year):
    '''
    计算西伯利亚高压
    '''
    print(np.shape(Field))
    lons = np.arange(0, 360, 2.5, dtype=float)
    lats = np.arange(90, -90 - 1, -2.5, dtype=float)

    lat1 = np.where(lats <= 20,True,False)
    lat2 = np.where(lats >= 15,True,False)
    lat = np.logical_and(lat1,lat2)

    lon1 = np.where(lons <= 100,True,False)
    lon2 = np.where(lons >= 80,True,False)
    lon = np.logical_and(lon1,lon2)
    print(lat.shape,lon.shape)

    Field2 = Field[:,:,lon]
    Field2 = Field2[:,lat,:]
    n1 = Field2.shape
    Field2 = Field2.reshape(n1[0],-1)
    Field2=np.mean(Field2,axis=1)
    Field2,a =dclim.mapstd(Field2)
    Field2  = Field2.flatten()
    print('c1=',Field2.shape)
    print('c2=',np.array(I_Year).shape)

    return Field2,np.array(I_Year)

def makeValueGridzWithMask(x, y, values, prj_path):
    from .parameter import get_param_value
    xi, yi = np.linspace(min(x), max(x), 200), np.linspace(min(y), max(y), 200)
    grid_x, grid_y = np.meshgrid(xi, yi)
    grid_z = scipy.interpolate.griddata((x, y), values, (grid_x, grid_y), method='linear')

    # structure = getParamValue('structure', prj_path)
    tunnel_thick = float(get_param_value('tc.tunnel.thick', prj_path))
    trap_thick = float(get_param_value('tc.trap.thick', prj_path))
    block_thick = float(get_param_value('tc.block.thick', prj_path))
    iso1_width = float(get_param_value('tc.iso1.width', prj_path))
    gate1_width = float(get_param_value('tc.gate1.width', prj_path))
    iso2_width = float(get_param_value('tc.iso2.width', prj_path))
    gate2_width = float(get_param_value('tc.gate2.width', prj_path))
    iso3_width = float(get_param_value('tc.iso3.width', prj_path))
    gate3_width = float(get_param_value('tc.gate3.width', prj_path))
    iso4_width = float(get_param_value('tc.iso4.width', prj_path))
    main_thick = tunnel_thick + trap_thick + block_thick

    mask_y = np.array(grid_y > main_thick)
    mask_x_gate1 = np.logical_and(grid_x > iso1_width, grid_x < iso1_width + gate1_width)
    mask_x_gate2 = np.logical_and(grid_x > iso1_width + gate1_width + iso2_width,
                                  grid_x < iso1_width + gate1_width + iso2_width + gate2_width)
    mask_x_gate3 = np.logical_and(grid_x > iso1_width + gate1_width + iso2_width + gate2_width + iso3_width,
                                  grid_x < iso1_width + gate1_width + iso2_width + gate2_width + iso3_width
                                  + gate3_width)
    mask_z = mask_y & (mask_x_gate1 | mask_x_gate2 | mask_x_gate3)
    grid_z_masked = np.ma.array(grid_z, mask=mask_z)
    return grid_z_masked

    def refine_Hessian(self, kpx, kpy, kps):
        """
        Refine the keypoint location based on a 3 point derivative, and delete
        non-coherent keypoints.

        :param kpx: x_pos of keypoint
        :param kpy: y_pos of keypoint
        :param kps: s_pos of keypoint
        :return: arrays of corrected coordinates of keypoints, values and
            locations of keypoints
        """
        curr = self.dogs[(kps, kpy, kpx)]
        nx = self.dogs[(kps, kpy, kpx + 1)]
        px = self.dogs[(kps, kpy, kpx - 1)]
        ny = self.dogs[(kps, kpy + 1, kpx)]
        py = self.dogs[(kps, kpy - 1, kpx)]
        ns = self.dogs[(kps + 1, kpy, kpx)]
        ps = self.dogs[(kps - 1, kpy, kpx)]

        nxny = self.dogs[(kps, kpy + 1, kpx + 1)]
        nxpy = self.dogs[(kps, kpy - 1, kpx + 1)]
        pxny = self.dogs[(kps, kpy + 1, kpx - 1)]
        pxpy = self.dogs[(kps, kpy - 1, kpx - 1)]

        nsny = self.dogs[(kps + 1, kpy + 1, kpx)]
        nspy = self.dogs[(kps + 1, kpy - 1, kpx)]
        psny = self.dogs[(kps - 1, kpy + 1, kpx)]
        pspy = self.dogs[(kps - 1, kpy - 1, kpx)]

        nxns = self.dogs[(kps + 1, kpy, kpx + 1)]
        nxps = self.dogs[(kps - 1, kpy, kpx + 1)]
        pxns = self.dogs[(kps + 1, kpy, kpx - 1)]
        pxps = self.dogs[(kps - 1, kpy, kpx - 1)]

        dx = (nx - px) / 2.0
        dy = (ny - py) / 2.0
        ds = (ns - ps) / 2.0
        dxx = (nx - 2.0 * curr + px)
        dyy = (ny - 2.0 * curr + py)
        dss = (ns - 2.0 * curr + ps)
        dxy = (nxny - nxpy - pxny + pxpy) / 4.0
        dxs = (nxns - nxps - pxns + pxps) / 4.0
        dsy = (nsny - nspy - psny + pspy) / 4.0
        det = -(dxs * dyy * dxs) + dsy * dxy * dxs + dxs * dsy * dxy - dss * dxy * dxy - dsy * dsy * dxx + dss * dyy * dxx
        K00 = dyy * dxx - dxy * dxy
        K01 = dxs * dxy - dsy * dxx
        K02 = dsy * dxy - dxs * dyy
        K10 = dxy * dxs - dsy * dxx
        K11 = dss * dxx - dxs * dxs
        K12 = dxs * dsy - dss * dxy
        K20 = dsy * dxy - dyy * dxs
        K21 = dsy * dxs - dss * dxy
        K22 = dss * dyy - dsy * dsy

        delta_s = -(ds * K00 + dy * K01 + dx * K02) / det
        delta_y = -(ds * K10 + dy * K11 + dx * K12) / det
        delta_x = -(ds * K20 + dy * K21 + dx * K22) / det
        peakval = curr + 0.5 * (delta_s * ds + delta_y * dy + delta_x * dx)
        mask = numpy.logical_and(numpy.logical_and(abs(delta_x) < self.tresh, abs(delta_y) < self.tresh), abs(delta_s) < self.tresh)
        return kpx + delta_x, kpy + delta_y, kps + delta_s, peakval, mask

 def _call_joint_genotypes(self, data, genotypes):
     normal = genotypes['normal']
     tumour = genotypes['tumour']
             
     normal_aa = (normal == 0)
     normal_ab = (normal == 1)
     normal_bb = (normal == 2)
         
     normal_var = np.logical_or(normal_ab, normal_bb)
     
     tumour_aa = (tumour == 0)
     tumour_ab = (tumour == 1)
     tumour_bb = (tumour == 2)
     
     tumour_var = np.logical_or(tumour_ab, tumour_bb)
     tumour_hom = np.logical_and(tumour_aa, tumour_bb)
     
     reference = np.logical_and(normal_aa, tumour_aa)
     germline = np.logical_and(normal_var, tumour_var)
     somatic = np.logical_and(normal_aa, tumour_var)
     loh = np.logical_and(normal_ab, tumour_hom)
     
     
     n = normal_aa.size
     joint_genotypes = 4 * np.ones((n,))
     
     joint_genotypes[reference] = 0
     joint_genotypes[germline] = 1
     joint_genotypes[somatic] = 2
     joint_genotypes[loh] = 3
     
     return joint_genotypes

def analyzeResult(x, accuracy, perturbAt=10000, movingAvg=True, smooth=True):
  if movingAvg:
    accuracy = movingAverage(accuracy, min(len(accuracy), 100))

  x = np.array(x)
  accuracy = np.array(accuracy)
  if smooth:
    # perform smoothing convolution
    mask = np.ones(shape=(100,))
    mask = mask/np.sum(mask)
    # extend accuracy vector to eliminate boundary effect of convolution
    accuracy = np.concatenate((accuracy, np.ones((200, ))*accuracy[-1]))
    accuracy = np.convolve(accuracy, mask, 'same')
    accuracy = accuracy[:len(x)]


  perturbAtX = np.where(x > perturbAt)[0][0]

  finalAccuracy = accuracy[perturbAtX-len(mask)/2]
  learnTime = min(np.where(np.logical_and(accuracy > finalAccuracy * 0.99,
                                          x < x[perturbAtX - len(mask)/2-1]))[0])
  learnTime = x[learnTime]

  finalAccuracyAfterPerturbation = accuracy[-1]
  learnTimeAfterPerturbation = min(np.where(
    np.logical_and(accuracy > finalAccuracyAfterPerturbation * 0.99,
                   x > x[perturbAtX + len(mask)]))[0])

  learnTimeAfterPerturbation = x[learnTimeAfterPerturbation] - perturbAt

  result = {"finalAccuracy": finalAccuracy,
            "learnTime": learnTime,
            "finalAccuracyAfterPerturbation": finalAccuracyAfterPerturbation,
            "learnTimeAfterPerturbation": learnTimeAfterPerturbation}
  return result

  def _compute_health_pill(self, x):
    x_clean = x[np.where(
        np.logical_and(
            np.logical_not(np.isnan(x)), np.logical_not(np.isinf(x))))]
    if np.size(x_clean):
      x_min = np.min(x_clean)
      x_max = np.max(x_clean)
      x_mean = np.mean(x_clean)
      x_var = np.var(x_clean)
    else:
      x_min = np.inf
      x_max = -np.inf
      x_mean = np.nan
      x_var = np.nan

    return np.array([
        1.0,  # Assume is initialized.
        np.size(x),
        np.sum(np.isnan(x)),
        np.sum(x == -np.inf),
        np.sum(np.logical_and(x < 0.0, x != -np.inf)),
        np.sum(x == 0.0),
        np.sum(np.logical_and(x > 0.0, x != np.inf)),
        np.sum(x == np.inf),
        x_min,
        x_max,
        x_mean,
        x_var,
        float(tf.as_dtype(x.dtype).as_datatype_enum),
        float(len(x.shape)),
    ] + list(x.shape))

def res2_true_and_false(hs, res, SV):
    'Organizes results into true positive and false positive sets'
    if not 'SV' in vars():
        SV = True
    #if not 'res' in vars():
        #res = qcx2_res[qcx]
    indx_samp = hs.indexed_sample_cx
    qcx = res.qcx
    cx2_score = res.cx2_score if SV else res.cx2_score
    unfilt_top_cx = np.argsort(cx2_score)[::-1]
    # Get top chip indexes and scores
    top_cx    = np.array(helpers.intersect_ordered(unfilt_top_cx, indx_samp))
    top_score = cx2_score[top_cx]
    # Get the true and false ground truth ranks
    qnx         = hs.tables.cx2_nx[qcx]
    if qnx <= 1:
        qnx = -1  # disallow uniden animals from being marked as true
    top_nx      = hs.tables.cx2_nx[top_cx]
    true_ranks  = np.where(np.logical_and(top_nx == qnx, top_cx != qcx))[0]
    false_ranks = np.where(np.logical_and(top_nx != qnx, top_cx != qcx))[0]
    # Construct the true positive tuple
    true_scores  = top_score[true_ranks]
    true_cxs     = top_cx[true_ranks]
    true_tup     = (true_cxs, true_scores, true_ranks)
    # Construct the false positive tuple
    false_scores = top_score[false_ranks]
    false_cxs    = top_cx[false_ranks]
    false_tup    = (false_cxs, false_scores, false_ranks)
    # Return tuples
    return true_tup, false_tup

def count_edges_within_band(a, b, band=3, rising=True):
    '''
    Counts the number of rising (or falling) edges match, within a sample band
    
    Params
    -------
    @param a, b: Arrays that will be compared
    @type a, b: Boolean array 
    @param band: The number of samples of tolerance
    @type band: float
    @param rising: Specify rising or falling edge
    @type rising: boolean 
    
    Returns
    -------
    @return: Count of matching edges, total true rising (or falling) edges
    @rtype: int
    '''
    if rising:
        a = np.r_[a[0], np.diff(a)]>0
        b = np.r_[b[0], np.diff(b)]>0
    else:
        a = np.r_[a[0], np.diff(a)]<0
        b = np.r_[b[0], np.diff(b)]<0

    total_edges = sum(a)
    result = np.logical_and(a, b)
    for offset in np.add(range(3),1):
        posoff = np.r_[[0]*offset, np.logical_and(a[:-offset], b[offset:])]
        negoff = np.r_[np.logical_and(a[offset:], b[:-offset]), [0]*offset]
        result = np.logical_or(result, posoff)
        result = np.logical_or(result, negoff)

    return sum(result), total_edges

def filter_params(io, rsh, rs, ee, isc):
    # Function filter_params identifies bad parameter sets. A bad set contains
    # Nan, non-positive or imaginary values for parameters; Rs > Rsh; or data
    # where effective irradiance Ee differs by more than 5% from a linear fit
    # to Isc vs. Ee

    badrsh = np.logical_or(rsh < 0., np.isnan(rsh))
    negrs = rs < 0.
    badrs = np.logical_or(rs > rsh, np.isnan(rs))
    imagrs = ~(np.isreal(rs))
    badio = np.logical_or(~(np.isreal(rs)), io <= 0)
    goodr = np.logical_and(~badrsh, ~imagrs)
    goodr = np.logical_and(goodr, ~negrs)
    goodr = np.logical_and(goodr, ~badrs)
    goodr = np.logical_and(goodr, ~badio)

    matrix = np.vstack((ee / 1000., np.zeros(len(ee)))).T
    eff = np.linalg.lstsq(matrix, isc)[0][0]
    pisc = eff * ee / 1000
    pisc_error = np.abs(pisc - isc) / isc
    # check for departure from linear relation between Isc and Ee
    badiph = pisc_error > .05

    u = np.logical_and(goodr, ~badiph)
    return u

def computeBinnedRates(
        spk_bins, times, step_size=5, bin_size=50, direction='forward'):

    if direction == 'forward':
        time_windows = np.c_[
            np.arange(times[0], times[-1] - bin_size + 1, step_size),
            np.arange(times[0], times[-1] - bin_size + 1,
                      step_size) + bin_size]
    elif direction == 'backward':
        time_windows = np.c_[
            np.arange(times[-1], times[0] + bin_size - 1,
                      -step_size) - bin_size,
            np.arange(times[-1], times[0] + bin_size - 1, -step_size)]

    n_trials = spk_bins.shape[0]
    n_windows = time_windows.shape[0]

    spk_counts = np.zeros([n_trials, n_windows])
    spk_counts[:] = np.nan

    for win_num in range(n_windows):
        if direction == 'forward':
            spk_counts[:, win_num] = np.sum(
                spk_bins[:, np.logical_and(
                    times >= time_windows[win_num][0],
                    times < time_windows[win_num][1])], axis=1)
            bin_centers = np.mean(time_windows, axis=1)
        elif direction == 'backward':
            spk_counts[:, -(win_num + 1)] = np.sum(
                spk_bins[:, np.logical_and(
                    times >= time_windows[win_num][0],
                    times < time_windows[win_num][1])], axis=1)
            bin_centers = np.mean(time_windows, axis=1)[::-1]
    spk_rates = spk_counts / bin_size * 1000
    return spk_rates, bin_centers

 def compare(self, result, t , noiseType = 'Z', sum = True, full=True):
     """Compares the discretization of this case with the one of an algorithm whose results are given in otherdisc. timeparam variable contains the variables for the discretization. Returns a dictionnary with the number of True positives, True negatives, False positives and False negatives"""
     #restrict comparation between Tb and Te
     
     otherdisc=result
     try:
         fulldisc=self.case[noiseType].discretize(t)
         intdisc=[int(b) for b in fulldisc]
         assert( len(otherdisc) == len(fulldisc) )
         #assert(not any([i==None for i in disc]))
     except AssertionError:
         print('something wrong in function of ', self)
     if full:
         disc=fulldisc
     else: 
         mask = np.logical_and(t >= self.case['Tb'], t <= self.case['Te'])
         otherdisc = result[mask]
         disc = fulldisc[mask]
         t=t[mask]
     retTF={}
     retTF['TP'] = np.logical_and(otherdisc,disc)
     retTF['TN'] = np.logical_and(np.logical_not(otherdisc), np.logical_not(disc))
     retTF['FP'] = np.logical_and(otherdisc, np.logical_not(disc))
     retTF['FN'] = np.logical_and(np.logical_not(otherdisc),  disc)
     if sum:
         for k, v in retTF.items():
             retTF[k]= int(v.sum())
     else:
         retTF['t'] = t
         retTF['disc'] = disc
     return retTF, intdisc

def plot_subplots(offsets, tag, n, window):
  fig, axarr = plt.subplots(4, sharex=True, **dict(figsize=(12,12)))
  sel1 = np.logical_and(offsets['mchirp'] > 0.0, offsets['mchirp'] <= 5.0)
  sel2 = np.logical_and(offsets['mchirp'] > 5.0, offsets['mchirp'] <= 10.0)
  sel3 = np.logical_and(offsets['mchirp'] > 10.0, offsets['mchirp'] <= 15.0)
  sel4 = np.logical_and(offsets['mchirp'] > 15.0, offsets['mchirp'] <= 100.0)
  
  sel = [sel1, sel2, sel3, sel4]
  labels = [r'$\mathcal{M}\ \leq\ 5M_{\odot}$', r'$5M_{\odot}\ <\ \mathcal{M}\ \leq\ 10M_{\odot}$'\
      ,r'$10M_{\odot}\ <\ \mathcal{M}\ \leq\ 15M_{\odot}$', r'$\mathcal{M}\ >\ 15M_{\odot}$']
  
  hist_min = np.min(offsets['offset'][sel[0]])
  hist_max = np.max(offsets['offset'][sel[0]])
  
  for i in xrange(len(axarr)):
    if not np.any(sel[i]):
      continue
    axarr[i].hist(offsets['offset'][sel[i]], histtype='step', bins=n,\
        range=[hist_min, hist_max], label=labels[i])
    axarr[i].set_yscale('log', nonposy='clip')
    axarr[i].set_xlim(-window/2, window/2)
    axarr[i].set_ylabel(r'N')
    axarr[i].grid(True)
    axarr[i].legend(loc="upper left", bbox_to_anchor=(1,1))

  axarr[3].set_xlabel(r'Offset  [Sec]')
  fig.subplots_adjust(hspace=0.5)
  plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)

  axarr[0].set_title( tag)
  plt.savefig(tag + '_subplotshistogram.png', bbox_inches='tight')
  #, bbox_extra_artists=(lgd,), bbox_inches='tight')
  plt.close()

    def close_gripper(self, lr, step_viewer=1, max_vel=.02, close_dist_thresh=0.004, grab_dist_thresh=0.005):
        print 'CLOSING GRIPPER'
        # generate gripper finger trajectory
        joint_ind = self.robot.GetJoint("%s_gripper_l_finger_joint" % lr).GetDOFIndex()
        start_val = self.robot.GetDOFValues([joint_ind])[0]
        print 'start_val: ', start_val
        # execute gripper finger trajectory
        dyn_bt_objs = [bt_obj for sim_obj in self.dyn_sim_objs for bt_obj in sim_obj.get_bullet_objects()]
        next_val = start_val
        while next_val:
            flr2finger_pts_grid = self._get_finger_pts_grid(lr)
            ray_froms, ray_tos = flr2finger_pts_grid['l'], flr2finger_pts_grid['r']

            # stop closing if any ray hits a dynamic object within a distance of close_dist_thresh from both sides
            next_vel = max_vel
            for bt_obj in dyn_bt_objs:
                from_to_ray_collisions = self.bt_env.RayTest(ray_froms, ray_tos, bt_obj)
                to_from_ray_collisions = self.bt_env.RayTest(ray_tos, ray_froms, bt_obj)
                rays_dists = np.inf * np.ones((len(ray_froms), 2))
                for rc in from_to_ray_collisions:
                    ray_id = np.argmin(np.apply_along_axis(np.linalg.norm, 1, ray_froms - rc.rayFrom))
                    rays_dists[ray_id, 0] = np.linalg.norm(rc.pt - rc.rayFrom)
                for rc in to_from_ray_collisions:
                    ray_id = np.argmin(np.apply_along_axis(np.linalg.norm, 1, ray_tos - rc.rayFrom))
                    rays_dists[ray_id, 1] = np.linalg.norm(rc.pt - rc.rayFrom)
                colliding_rays_inds = np.logical_and(rays_dists[:, 0] != np.inf, rays_dists[:, 1] != np.inf)
                if np.any(colliding_rays_inds):
                    rays_dists = rays_dists[colliding_rays_inds, :]
                    if np.any(np.logical_and(rays_dists[:, 0] < close_dist_thresh,
                                             rays_dists[:, 1] < close_dist_thresh)):
                        next_vel = 0
                    else:
                        next_vel = np.minimum(next_vel, np.min(rays_dists.sum(axis=1)))
            if next_vel == 0:
                break
            next_val = np.maximum(next_val - next_vel, 0)

            self.robot.SetDOFValues([next_val], [joint_ind])
            self.step()
            if self.viewer and step_viewer:
                self.viewer.Step()
        handles = []
        # add constraints at the points where a ray hits a dynamic link within a distance of grab_dist_thresh
        for bt_obj in dyn_bt_objs:
            from_to_ray_collisions = self.bt_env.RayTest(ray_froms, ray_tos, bt_obj)
            to_from_ray_collisions = self.bt_env.RayTest(ray_tos, ray_froms, bt_obj)
            
            for i in range(ray_froms.shape[0]):
                self.viewer.Step()
            ray_collisions = [rc for rcs in [from_to_ray_collisions, to_from_ray_collisions] for rc in rcs]

            for rc in ray_collisions:
                if rc.link == bt_obj.GetKinBody().GetLink('rope_59'):
                    self.viewer.Step()
                if np.linalg.norm(rc.pt - rc.rayFrom) < grab_dist_thresh:
                    link_tf = rc.link.GetTransform()
                    link_tf[:3, 3] = rc.pt
                    self._add_constraints(lr, rc.link, link_tf)
        if self.viewer and step_viewer:
            self.viewer.Step()