def occupancy_indicator(ch2bins, peak_sets):
    """
This function determines the occupancy indicator for each pair of bin and peak
set.

ch2bins = { chrom: [(start, end, summit, ID), ...], ... }, where bins for each
chromatin have been sorted based on (start, end), and ID's must range from 1
through the total number of bins. Each start coordinate must be strictly less
than the corresponding end site. Summit is not used in this function.

peak_sets is a list of peak sets, each having the format:
{ chrom: [(start, end, ...), ...], ... }, where peaks for each chromatin have
been sorted based on (start, end), and each start coordinate must be strictly
less than the corresponding end site.

Note: for ch2bins or each single peak set of peak_sets, there mustn't be a
region that is enclosed by another one.

Return value = [[0, 1, ...], ...]. For a bin of ID i, ret[i - 1][j] = 1 iff the
bin is occupied by the peak set at index j, which means the bin's middle point
falls within some peak of the set.

"""
    bin_num = sum(map(len, ch2bins.values()))
    sample_num = len(peak_sets)
    occupy = [None] * bin_num
    for i in range(bin_num):
        occupy[i] = [0] * sample_num
    for j, peak_set in enumerate(peak_sets):
        for ch in ch2bins:
            if ch not in peak_set:
                continue
            bins = ch2bins[ch]
            peaks = peak_set[ch]
            pind = 0
            for start, end, temp, id in bins:
                mid = (start + end - 1) // 2
                while pind < len(peaks):
                    s, e = peaks[pind][:2]
                    if mid < s:
                        break
                    if mid < e:
                        occupy[id - 1][j] = 1
                        break
                    pind += 1
                else:
                    break
    return occupy
