Implementing the DWT
Treating boundaries correctly is non trivial for wavelets and filter banks in general. Symmetric/anti-symmetric filters (such as with biorthogonal wavelets as the 5/3 or the 9/7) help a little. Yet, with dyadic sub-sampling, without playing too heavily on polyphase matrices and symmetries, I would suggest the following procedure:


choose a level $L$ for the iterated decomposition 
find an image size that can be divided by $2^L$
extend the original image to that size (symmetry/anti-symmetry, or zero-padding but not the best)
extend the smallest filter with zero-padding to match the longest filter size, so both filters have the same padded length.


Cleverer/sparser exist, but they are more involved. Now the number of coefficients ought to be similar. As for the Python, this is more "language-related", and not adapted to SE.DSP. IN 1D, the lifting scheme can improve the computation burden by a factor of 2 in terms of basioc operations.