I think what you want to do can be accomplished using bedtools map.

You could either use bedops --partition followed by bedtools map, or simply use bedtools unionbedg and post-process its output. If you're using pybedtools, the latter option is probably easiest. There may also be a way to do all of this with bedops.

I think you can do this with a one-liner with BEDOPS tools:

$ bedops --partition elements.bed | bedmap --echo --mean --delim '\t' - elements.bed > answer.bedgraph

The elements.bed file should follow UCSC conventions and have its score data in the fifth column.

The way this works is that elements in elements.bed are partitioned into disjoint ranges, and then the disjoint ranges (the "-" in bedmap) are mapped back against the elements, with arithmetic means calculated from the score values of overlapping elements. The delim operand should put results into bedgraph format.

The only requirement for BEDOPS is that inputs are sorted prep the file, if necessary:

$ sort-bed elements.unsorted.bed > elements.bed

Output from BEDOPS tools is in correct sort order; you only need to sort once. Sorting with sort-bed is generally faster than alternatives and is non-lossy.

If you need to do all this in Python, you could use subprocess.

Of course, this requires the willingness to use BEDOPS, which not all want to do for reasons. But I'd suggest using whatever solves the problem efficiently. At the end of the day, an external dependency is not much different from an internal dependency - if it doesn't come with the language or the shell, it's a third-party API that needs installation and maintenance, either way.