dtrack_utils_testing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "from libc.math cimport abs, sqrt, ceil, floor, log, log2, acos, cos\n",
    "from numpy cimport ndarray\n",
    "import numpy as np\n",
    "from numpy import ones, zeros, int32, float32, uint8, fromstring\n",
    "from numpy import sort, empty, array, arange, concatenate, searchsorted\n",
    "from numpy import minimum, maximum,divide, std, mean\n",
    "from numpy import min as nmin\n",
    "from numpy import max as nmax\n",
    "import h5py as h\n",
    "\n",
    "def pairRegionsIntersection(ndarray[int, ndim=2] pairs,\n",
    "                            ndarray[int, ndim=2] regions,\n",
    "                            exclude=False, allow_partial=False,\n",
    "                            bool_out=False):\n",
    "    '''\n",
    "    Given an list of pairs of the form [a,b] where a<b and regions of the form [c,d] we want to return the indices\n",
    "    (within the pair list) of pairs which overlap with the regions. This will be used to exclude certain regions of the\n",
    "    genome if needed.\n",
    "     \n",
    "    Arguments:\n",
    "    \n",
    "    - pairs: (N,2) shape array where each row details a pair.\n",
    "    - regions: (M,2) shape array where each row details a region.\n",
    "    - exclude: Boolean. If true, then return the indices of pairs which don't overlap with any of the regions.\n",
    "    - allow_partial: Boolean. If true then include pair indices which only partially overlap with a given region.\n",
    "                     Here we consider the 'partial' to apply to the pair. That is, if allow_partial == False then\n",
    "                     we must satisfy either:\n",
    "                         - the entire pair is included in a region if exclude == False\n",
    "                         - The entire pair must not be included in a region if exclude == True\n",
    "                     If allow_partial == True then we must satisfy:\n",
    "                         - At least part of the pair touches at least one region if exclude == False\n",
    "                         - At least part of the pair touches no regions if exclude == True\n",
    "                      \n",
    "    Returns:\n",
    "    \n",
    "    - indices: The indices of the pairs in the pair array which satisfy the overlapping conditions.\n",
    "    '''\n",
    "    cdef int i, j, k, a, b\n",
    "    cdef int exc = int(exclude)\n",
    "    cdef int partial = int(allow_partial)\n",
    "    cdef int bool_o = int(bool_out)\n",
    "    cdef int ni = 0\n",
    "    cdef int np = len(pairs)\n",
    "    cdef int nr = len(regions)\n",
    "    cdef ndarray[int, ndim=1] indices = zeros(np, int32)\n",
    "    cdef ndarray[int, ndim=1] indices_out = empty(np, int32)\n",
    "    cdef ndarray[int, ndim=1] order = array(regions[:,0].argsort(), int32)\n",
    "    cdef ndarray[int, ndim=2] regs = empty((nr+1,2), int32)\n",
    "    cdef minpoint = minimum(nmin(pairs), nmin(regions))\n",
    "    cdef maxpoint = maximum(nmax(pairs), nmax(regions))\n",
    "\n",
    "    if exc:\n",
    "        for i in range(nr):\n",
    "            regs[i,1] = regions[order[i],0]\n",
    "        regs[0,0] = minpoint - 1\n",
    "        regs[-1,1] = maxpoint + 1\n",
    "    else:\n",
    "        regs = regions\n",
    "    \n",
    "    \n",
    "    for i in range(np):\n",
    "        for k in range(nr):\n",
    "            j = order[k]\n",
    "            if (regs[j,1] > pairs[i,0]) and (pairs[i,1] > regs[j,0]):\n",
    "                #Pair is at least partially overlapping with the region\n",
    "                if partial:\n",
    "                    indices[i] = 1\n",
    "                    if bool_o:\n",
    "                        return 1\n",
    "                    continue\n",
    "                elif (regs[j,1] >= pairs[i,1]) and (pairs[i,0] >= regs[j,0]):\n",
    "                    #Pair is entirely containing within the region\n",
    "                    indices[i] = 1\n",
    "                    if bool_o:\n",
    "                        return 1\n",
    "                    continue\n",
    "    \n",
    "    if bool_o:\n",
    "        return 0\n",
    "    \n",
    "    for i in range(np):\n",
    "        if indices[i] == 1:\n",
    "            indices_out[ni] = i\n",
    "            ni +=1\n",
    "    \n",
    "    return indices_out[:ni]\n",
    "\n",
    "def multi_pairRegionsIntersection(ndarray[int, ndim=3] pairs,\n",
    "                                  ndarray[int, ndim=3] regions,\n",
    "                                  exclude=False, allow_partial=False, indices = False):\n",
    "    '''\n",
    "    Given an list of pairs of the form [a,b] where a<b and regions of the form [c,d] we want to work pairwise to\n",
    "    return the region indices of the regions which overlap with each pair.\n",
    "     \n",
    "    Arguments:\n",
    "    \n",
    "    - pairs: (N,2) shape array where each row details a pair.\n",
    "    - regions: (M,2) shape array where each row details a region.\n",
    "    - exclude: Boolean. If true, then return the indices of pairs which don't overlap with any of the regions.\n",
    "    - allow_partial: Boolean. If true then include pair indices which only partially overlap with a given region.\n",
    "                     Here we consider the 'partial' to apply to the pair. That is, if allow_partial == False then\n",
    "                     we must satisfy either:\n",
    "                         - the entire pair is included in a region if exclude == False\n",
    "                         - The entire pair must not be included in a region if exclude == True\n",
    "                     If allow_partial == True then we must satisfy:\n",
    "                         - At least part of the pair touches at least one region if exclude == False\n",
    "                         - At least part of the pair touches no regions if exclude == True\n",
    "                      \n",
    "    Returns:\n",
    "    \n",
    "    - out: A maximum of an (N*M,2) shape array where each column details a pair group index and a regions group\n",
    "           index where than pair group and that region group overlap. This is essentially gonna be in COO\n",
    "           format\n",
    "    '''\n",
    "    cdef int i, j, k, np, nr,nl\n",
    "    cdef int overlap\n",
    "    cdef int ind = int(indices)\n",
    "    cdef int maxnp = len(pairs[0,:,0]) \n",
    "    cdef int maxnr = len(regions[0,:,0])\n",
    "    cdef int minpairs = nmin(pairs)\n",
    "    cdef int minregions = nmin(regions)\n",
    "    cdef int npairgroups = len(pairs)\n",
    "    cdef int nreggroups = len(regions)\n",
    "    cdef ndarray[int, ndim=2] out = zeros((npairgroups*nreggroups,2), int32)\n",
    "    \n",
    "    for i in range(npairgroups):\n",
    "        print(\"#################################\")\n",
    "        print(\"pair {}\".format(i))\n",
    "        print(\"Unedited pairs:\")\n",
    "        print(pairs[i,:,:])\n",
    "        np = maxnp\n",
    "        for k in range(maxnp):\n",
    "            if pairs[i,k,0] == minpairs:\n",
    "                np = k\n",
    "                break\n",
    "        for j in range(nreggroups):\n",
    "            print(\"#################################\")\n",
    "            print(\"region {}\".format(j))\n",
    "            print(\"Unedited regions:\")\n",
    "            print(regions[j,:,:])\n",
    "            nr = maxnr\n",
    "            for k in range(maxnr):\n",
    "                if regions[j,k,0] == minregions:\n",
    "                    nr = k\n",
    "                    break\n",
    "            print(\"******************************\")\n",
    "            print(\"Edited pairs:\")\n",
    "            print(pairs[i,:np,:])\n",
    "            print(\"Edited regions:\")\n",
    "            print(regions[j,:nr,:])\n",
    "            overlap = pairRegionsIntersection(pairs[i,:np,:],\n",
    "                                              regions[j,:nr,:],\n",
    "                                              allow_partial = allow_partial,\n",
    "                                              exclude = exclude,\n",
    "                                              bool_out = True)\n",
    "            print(\"Output: {}\".format(overlap))\n",
    "            if overlap == 1:\n",
    "                out[nl,0] = i\n",
    "                out[nl,1] = j\n",
    "                nl +=1\n",
    "    \n",
    "    return out[:nl,:]\n",
    "                \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "pairs = np.array([[[10,50],\n",
    "                   [51,52],\n",
    "                   [55,65],\n",
    "                   [70,100]],[[5,20],\n",
    "                              [21,22],\n",
    "                              [-1,-1],\n",
    "                              [-1,-1]],[[110,130],\n",
    "                                        [-1,-1],\n",
    "                                        [-1,-1],\n",
    "                                        [-1,-1]]]).astype('int32')\n",
    "regions = np.array([[[0,5],\n",
    "                     [7,9],\n",
    "                     [-1,-1],\n",
    "                     [-1,-1],\n",
    "                     [-1,-1]],[[2,7],\n",
    "                               [10,22],\n",
    "                               [7,17],\n",
    "                               [21,100],\n",
    "                               [105,151]],[[10,20],\n",
    "                                           [21,22],\n",
    "                                           [34,47],\n",
    "                                           [51,100],\n",
    "                                           [101,130]],[[10,20],\n",
    "                                                       [30,40],\n",
    "                                                       [50,60],\n",
    "                                                       [-1,-1],\n",
    "                                                       [-1,-1]]]).astype('int32')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 4, 2)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#################################\n",
      "pair 0\n",
      "Unedited pairs:\n",
      "[[ 10  50]\n",
      " [ 51  52]\n",
      " [ 55  65]\n",
      " [ 70 100]]\n",
      "#################################\n",
      "region 0\n",
      "Unedited regions:\n",
      "[[ 0  5]\n",
      " [ 7  9]\n",
      " [-1 -1]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 10  50]\n",
      " [ 51  52]\n",
      " [ 55  65]\n",
      " [ 70 100]]\n",
      "Edited regions:\n",
      "[[0 5]\n",
      " [7 9]]\n",
      "Output: 0\n",
      "#################################\n",
      "region 1\n",
      "Unedited regions:\n",
      "[[  2   7]\n",
      " [ 10  22]\n",
      " [  7  17]\n",
      " [ 21 100]\n",
      " [105 151]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 10  50]\n",
      " [ 51  52]\n",
      " [ 55  65]\n",
      " [ 70 100]]\n",
      "Edited regions:\n",
      "[[  2   7]\n",
      " [ 10  22]\n",
      " [  7  17]\n",
      " [ 21 100]\n",
      " [105 151]]\n",
      "Output: 1\n",
      "#################################\n",
      "region 2\n",
      "Unedited regions:\n",
      "[[ 10  20]\n",
      " [ 21  22]\n",
      " [ 34  47]\n",
      " [ 51 100]\n",
      " [101 130]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 10  50]\n",
      " [ 51  52]\n",
      " [ 55  65]\n",
      " [ 70 100]]\n",
      "Edited regions:\n",
      "[[ 10  20]\n",
      " [ 21  22]\n",
      " [ 34  47]\n",
      " [ 51 100]\n",
      " [101 130]]\n",
      "Output: 1\n",
      "#################################\n",
      "region 3\n",
      "Unedited regions:\n",
      "[[10 20]\n",
      " [30 40]\n",
      " [50 60]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 10  50]\n",
      " [ 51  52]\n",
      " [ 55  65]\n",
      " [ 70 100]]\n",
      "Edited regions:\n",
      "[[10 20]\n",
      " [30 40]\n",
      " [50 60]]\n",
      "Output: 1\n",
      "#################################\n",
      "pair 1\n",
      "Unedited pairs:\n",
      "[[ 5 20]\n",
      " [21 22]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "#################################\n",
      "region 0\n",
      "Unedited regions:\n",
      "[[ 0  5]\n",
      " [ 7  9]\n",
      " [-1 -1]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 5 20]\n",
      " [21 22]]\n",
      "Edited regions:\n",
      "[[0 5]\n",
      " [7 9]]\n",
      "Output: 0\n",
      "#################################\n",
      "region 1\n",
      "Unedited regions:\n",
      "[[  2   7]\n",
      " [ 10  22]\n",
      " [  7  17]\n",
      " [ 21 100]\n",
      " [105 151]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 5 20]\n",
      " [21 22]]\n",
      "Edited regions:\n",
      "[[  2   7]\n",
      " [ 10  22]\n",
      " [  7  17]\n",
      " [ 21 100]\n",
      " [105 151]]\n",
      "Output: 1\n",
      "#################################\n",
      "region 2\n",
      "Unedited regions:\n",
      "[[ 10  20]\n",
      " [ 21  22]\n",
      " [ 34  47]\n",
      " [ 51 100]\n",
      " [101 130]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 5 20]\n",
      " [21 22]]\n",
      "Edited regions:\n",
      "[[ 10  20]\n",
      " [ 21  22]\n",
      " [ 34  47]\n",
      " [ 51 100]\n",
      " [101 130]]\n",
      "Output: 1\n",
      "#################################\n",
      "region 3\n",
      "Unedited regions:\n",
      "[[10 20]\n",
      " [30 40]\n",
      " [50 60]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[ 5 20]\n",
      " [21 22]]\n",
      "Edited regions:\n",
      "[[10 20]\n",
      " [30 40]\n",
      " [50 60]]\n",
      "Output: 0\n",
      "#################################\n",
      "pair 2\n",
      "Unedited pairs:\n",
      "[[110 130]\n",
      " [ -1  -1]\n",
      " [ -1  -1]\n",
      " [ -1  -1]]\n",
      "#################################\n",
      "region 0\n",
      "Unedited regions:\n",
      "[[ 0  5]\n",
      " [ 7  9]\n",
      " [-1 -1]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[110 130]]\n",
      "Edited regions:\n",
      "[[0 5]\n",
      " [7 9]]\n",
      "Output: 0\n",
      "#################################\n",
      "region 1\n",
      "Unedited regions:\n",
      "[[  2   7]\n",
      " [ 10  22]\n",
      " [  7  17]\n",
      " [ 21 100]\n",
      " [105 151]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[110 130]]\n",
      "Edited regions:\n",
      "[[  2   7]\n",
      " [ 10  22]\n",
      " [  7  17]\n",
      " [ 21 100]\n",
      " [105 151]]\n",
      "Output: 1\n",
      "#################################\n",
      "region 2\n",
      "Unedited regions:\n",
      "[[ 10  20]\n",
      " [ 21  22]\n",
      " [ 34  47]\n",
      " [ 51 100]\n",
      " [101 130]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[110 130]]\n",
      "Edited regions:\n",
      "[[ 10  20]\n",
      " [ 21  22]\n",
      " [ 34  47]\n",
      " [ 51 100]\n",
      " [101 130]]\n",
      "Output: 1\n",
      "#################################\n",
      "region 3\n",
      "Unedited regions:\n",
      "[[10 20]\n",
      " [30 40]\n",
      " [50 60]\n",
      " [-1 -1]\n",
      " [-1 -1]]\n",
      "******************************\n",
      "Edited pairs:\n",
      "[[110 130]]\n",
      "Edited regions:\n",
      "[[10 20]\n",
      " [30 40]\n",
      " [50 60]]\n",
      "Output: 0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0, 1],\n",
       "       [0, 2],\n",
       "       [0, 3],\n",
       "       [1, 1],\n",
       "       [1, 2],\n",
       "       [2, 1],\n",
       "       [2, 2]], dtype=int32)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multi_pairRegionsIntersection(pairs, regions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "from libc.math cimport abs, sqrt, ceil, floor, log, log2, acos, cos\n",
    "from numpy cimport ndarray\n",
    "import numpy as np\n",
    "from numpy import ones, zeros, int32, float32, uint8, fromstring\n",
    "from numpy import sort, empty, array, arange, concatenate, searchsorted\n",
    "from numpy import minimum, maximum,divide, std, mean\n",
    "from numpy import min as nmin\n",
    "from numpy import max as nmax\n",
    "import h5py as h\n",
    "\n",
    "def non_overlapping(ndarray[int, ndim=2] regions):\n",
    "    '''\n",
    "    Given a list of regions which may be overlapping we return a sorted list of regions which aren't overlapping\n",
    "    arguments:\n",
    "    - regions: An (N,2) shape array with each row detailing the start and end of a region\n",
    "    \n",
    "    outputs:\n",
    "    - nonoverlapping: An (M,2) shape array where each row details a covered region. The regions in the output array\n",
    "                      will be sorted by start position and non-overlapping with each other. \n",
    "    '''\n",
    "    cdef int idx, jdx, nr = len(regions)\n",
    "    cdef ndarray[int, ndim=2] outregs = zeros((nr,2),int32)\n",
    "    cdef ndarray[int, ndim=2] values = zeros((2*regions.shape[0],2), int32)\n",
    "    for idx in range(nr):\n",
    "        values[idx,0] = regions[idx,0]\n",
    "        values[idx,1] = 1\n",
    "        values[idx+nr,0] = regions[idx,1]\n",
    "        values[idx+nr,1] = -1\n",
    "    cdef ndarray[int, ndim=1] order = array(values[:,0].argsort(), int32)\n",
    "    \n",
    "    cdef int total = 0\n",
    "    cdef int new_point = 0\n",
    "    cdef int valtype = 0\n",
    "    cdef int current_start = 0\n",
    "    cdef int current_end = 0\n",
    "    cdef int out_idx = 0\n",
    "    for idx in np.arange(2*nr):\n",
    "        jdx = order[idx]\n",
    "        new_point, valtype = values[jdx]\n",
    "        total += valtype\n",
    "        if valtype <0:\n",
    "            current_end = new_point\n",
    "\n",
    "        if (total == 1) and valtype > 0:\n",
    "            current_start = new_point\n",
    "                \n",
    "        if total == 0:\n",
    "            outregs[out_idx,0] = current_start\n",
    "            outregs[out_idx,1] = current_end\n",
    "            out_idx += 1\n",
    "            \n",
    "    return outregs[:out_idx,:]\n",
    "                \n",
    "\n",
    "\n",
    "def pairRegionsOverlap(ndarray[int, ndim=2] pairs,\n",
    "                       ndarray[int, ndim=2] regions,\n",
    "                       exclude=False):\n",
    "\n",
    "    '''\n",
    "    Given some (N,2) shape array of pairs a,b with a<b and an (M,2) shape array of regions c,d, with c<d, return a new\n",
    "    (R,2) shape array of overlaps e,f with e<f such that each overlap is the intersection of the sets P and R where P\n",
    "    is the union of the pairs and R is the union of the regions. \n",
    "    \n",
    "    Arguments:\n",
    "    \n",
    "    - pairs: (N,2) shape array where each row details a pair.\n",
    "    - regions: (M,2) shape array where each row details a region.\n",
    "    - exclude: Boolean. If True then the function returns a new set of overlaps composed of the overlaps in the set\n",
    "               P union ~R. \n",
    "               \n",
    "    Returns:\n",
    "    - overlap_bounds: An (R,2) array detailing the set of intervals making up the overlaps between the pairs and\n",
    "                      regions\n",
    "    '''\n",
    "    \n",
    "    pairs = non_overlapping(pairs)\n",
    "    regions = non_overlapping(regions)\n",
    "    cdef int i, j, k\n",
    "    cdef int exc = int(exclude)\n",
    "    #totaloverlaps index counter\n",
    "    cdef int tnio = 0\n",
    "    #number of pairs\n",
    "    cdef int np = len(pairs)\n",
    "    #number of regions\n",
    "    cdef int nr = len(regions)\n",
    "    #order array for the regions\n",
    "    cdef ndarray[int, ndim=1] order = array(regions[:,0].argsort(), int32)\n",
    "    #Output regions array (maximum number of possible regions would be np + nr -1\n",
    "    cdef ndarray[int, ndim=2] overlap_bounds = empty(((np + nr),2), int32)\n",
    "    cdef ndarray[int, ndim=1] current_pair = empty(2, int32)\n",
    "    \n",
    "    for i in range(np):\n",
    "        current_pair = pairs[i,:]\n",
    "        \n",
    "        if current_pair[1] < regions[order[0],0]:\n",
    "            if exc:\n",
    "                overlap_bounds[tnio,:] = current_pair\n",
    "                tnio += 1\n",
    "                \n",
    "            continue\n",
    "        \n",
    "        if current_pair[0] > regions[order[nr-1],1]:\n",
    "            if exc:\n",
    "                overlap_bounds[tnio,:] = current_pair\n",
    "                tnio += 1\n",
    "                \n",
    "            continue\n",
    "            \n",
    "        for k in range(nr):\n",
    "            j = order[k]\n",
    "            #print(i,current_pair, j, regions[j], tnio)\n",
    "            if regions[j,0] < current_pair[1] and regions[j,1] > current_pair[0]:\n",
    "                #There is an overlap between this region and this pair chunk\n",
    "                if regions[j,0] <= current_pair[0] and current_pair[0] < regions[j,1]:\n",
    "                    #Start of the current pair chunk is in the region\n",
    "                    if regions[j,1] >= current_pair[1]:\n",
    "                        #Pair chunk entirely contained with a region \n",
    "                        if exc:\n",
    "                            #If exclude then we know that this\n",
    "                            #pair can't be included - it is contained within a \n",
    "                            #region so break the current loop and move onto the \n",
    "                            #next pair\n",
    "                            break\n",
    "                        else:\n",
    "                            #If include then we want to add this pair chunk to the \n",
    "                            #list of output pair chunks\n",
    "                            overlap_bounds[tnio] = current_pair\n",
    "                            tnio += 1\n",
    "                            continue\n",
    "                    else:\n",
    "                        #Start of pair chunk is in the region, end isnt\n",
    "                        if exc:\n",
    "                            #If exclude then we just shave off the bit of this\n",
    "                            #pair that touches the reegion and continue to the \n",
    "                            #next region\n",
    "                            current_pair[0] = regions[j,1]\n",
    "                            continue\n",
    "                        else:\n",
    "                            #If include then we add this overlap to the \n",
    "                            #output and shave off the bit of the pair that overlaps\n",
    "                            #with the region and then continue on to the next region\n",
    "                            overlap_bounds[tnio,0] = current_pair[0]\n",
    "                            overlap_bounds[tnio,1] = regions[j,1]\n",
    "                            \n",
    "                            tnio += 1\n",
    "                            \n",
    "                            current_pair[0] = regions[j,1]\n",
    "                            continue\n",
    "                if regions[j,0] > current_pair[0]:\n",
    "                    #Start of the current pair chunk is lower than the start of\n",
    "                    #the region\n",
    "                    if regions[j,1] < current_pair[1]:\n",
    "                        #End of the current pair chunk is higher than the start of the\n",
    "                        #regions - i.e. region entirely contained within current\n",
    "                        #pair chunk\n",
    "                        if exc:\n",
    "                            #Since the regions are sorted we know that the lower portion\n",
    "                            #of the pair which is cutoff can't overlap with any other \n",
    "                            #region so we add it to the output\n",
    "                            overlap_bounds[tnio,0] = current_pair[0]\n",
    "                            overlap_bounds[tnio,1] = regions[j,0]\n",
    "                            tnio += 1\n",
    "                            \n",
    "                            #We then shave the current pair to only include the bit\n",
    "                            #which didn't overlap with the current region above\n",
    "                            current_pair[0] = regions[j,1]\n",
    "                            continue\n",
    "                        else:\n",
    "                            #We include the overlapping region\n",
    "                            overlap_bounds[tnio,:] = regions[j,:]\n",
    "                            tnio += 1\n",
    "                            \n",
    "                            #Then shave the current pair\n",
    "                            current_pair[0] = regions[j,1]\n",
    "                            continue\n",
    "                    else:\n",
    "                        #Current pair end must overlap with the region\n",
    "                        if exc:\n",
    "                            overlap_bounds[tnio,0] = current_pair[0]\n",
    "                            overlap_bounds[tnio,1] = regions[j,0]\n",
    "                            tnio +=1\n",
    "                            #Since the regions are sorted we know that\n",
    "                            #we've already got all the info out\n",
    "                            #of this pair\n",
    "                            current_pair[0] = regions[j,0]\n",
    "                            break\n",
    "                        else:\n",
    "                            overlap_bounds[tnio,0] = regions[j,0]\n",
    "                            overlap_bounds[tnio,1] = current_pair[1]\n",
    "                            tnio += 1\n",
    "                            #Since the regions are sorted we know that\n",
    "                            #there won't be any more of this pair that\n",
    "                            #is overlapping any regions so continue\n",
    "                            #to the next pair\n",
    "                            current_pair[0] = regions[j,0]\n",
    "                            break\n",
    "            elif exc and regions[j,0] > current_pair[1]:\n",
    "                overlap_bounds[tnio,:] = current_pair \n",
    "                tnio += 1\n",
    "                break\n",
    "                    \n",
    "    \n",
    "    return overlap_bounds[:tnio]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs = np.array([[10,20],\n",
    "                   [30,40],\n",
    "                   [50,60],\n",
    "                   [70,90]]).astype('int32')\n",
    "regions = np.array([[15,25],\n",
    "                    [26,35],\n",
    "                    [38,55],\n",
    "                    [55,57],\n",
    "                    [75,151]]).astype('int32')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 15],\n",
       "       [35, 38],\n",
       "       [57, 60],\n",
       "       [70, 75]], dtype=int32)"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairRegionsOverlap(pairs,regions, exclude = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_root = \"/disk1/dh486/Data\"\n",
    "chip_paths = {'H3K4me3':'{}/ChIPseq/0hrs/peaks_bed/old_haploid/GSEXXXXXLaue_H3K4me3_hapmESC_peaks.npz'.format(data_root),\n",
    "              'H3K27me3':'{}/ChIPseq/0hrs/peaks_bed/old_haploid/GSEXXXXXLaue_H3K27me3_hapmESC_peaks.npz'.format(data_root),\n",
    "              'H3K27ac':'{}/ChIPseq/0hrs/peaks_bed/GSE56098_H3K27ac_ESC_peaks.npz'.format(data_root),\n",
    "              'H3K36me3':'{}/ChIPseq/0hrs/peaks_bed/old_haploid/GSEXXXXXLaue_H3K36me3_hapmESC_peaks.npz'.format(data_root)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = dict(np.load(chip_paths['H3K4me3']))['dataTrack/regions/H3K4me3/1']\n",
    "x2 = dict(np.load(chip_paths['H3K27me3']))['dataTrack/regions/H3K27me3/1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  3670867,   3671078],\n",
       "       [  3671153,   3671499],\n",
       "       [  3671734,   3671914],\n",
       "       ...,\n",
       "       [195176428, 195176645],\n",
       "       [195213000, 195213215],\n",
       "       [195241698, 195241961]], dtype=int32)"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairRegionsOverlap(x1.astype('int32'),x2.astype('int32'), exclude = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}