code: def subtract_ranges(r1s,r2s,already_sorted=False):
    from seqtools.stream import MultiLocusStream
    if not already_sorted:
      r1s = merge_ranges(r1s)
      r2s = merge_ranges(r2s)
    outputs = []
    mls = MultiLocusStream([BedArrayStream(r1s),BedArrayStream(r2s)])
    tot1 = 0
    tot2 = 0
    for loc in mls:
      v = loc.payload
      [beds1,beds2] =v
      beds1 = beds1[:]
      beds2 = beds2[:]
      if len(beds1)==0:
        continue
      if len(beds2)==0:
        outputs += beds1
        continue
      mapping = {}
      for i in range(0,len(beds1)):
        mapping[i] = []
      beds2min = 0
      beds2max = len(beds2)
      for i in range(0,len(beds1)):
        for j in range(beds2min,beds2max):
          cmpval = beds1[i].cmp(beds2[j])
          if cmpval == -1:
            beds2min = j+1
          elif cmpval == 0:
            mapping[i].append(j)
          else:
            break
      for i in range(0,len(beds1)):
        if len(mapping[i])==0: outputs += beds1
        else:
          outputs += subtract_range_array(beds1[i],[beds2[j] for j in mapping[i]],is_sorted=True)
    return merge_ranges(outputs)