Plotting randomly stacked cubes in 3D- mplot3d?

Question

I am exploring random stackings of cubes.

I started with 2D and can generate random packings within a rectangle like this: Now I have the code to generalize the stacking to 3D, but I am struggling to generalize the visualization. An example data set is, filling a 3x3x3 cube with 1x1x1 and 2x2x2 cubes,

#the coordinates of a corner vertex of the 19 1x1x1 cubes
x1 = [1, 0, 2, 0, 0, 0, 2, 1, 0, 1, 2, 2, 0, 0, 0, 2, 0, 1, 1]
y1 = [1, 1, 0, 2, 0, 0, 2, 2, 2, 0, 1, 0, 1, 2, 1, 0, 0, 0, 0]
z1 = [2, 1, 1, 0, 1, 2, 2, 2, 2, 1, 2, 0, 0, 1, 2, 2, 0, 0, 2]
#the coordinates of a corner vertex of the 1 2x2x2 cube 
x2 = [1]
y2 = [1]
z2 = [0]
# I believe the random filling is working because
# the total volumes equal:  19 + 2**3 = 3**3

#I would like to start with the lists  
X = [x1,x2]
Y = [y1,y2]
Z = [z1,z2]
sizes = [1,2]
#because I want to generalize the visualization to n sizes

So far, all I have the knowledge to do is plot a 3D scatter of the data

from matplotlib import pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for ii in range(len(sizes)): 
    ax.scatter(X[ii],Y[ii],Z[ii])
plt.show()

I would like to make a plot more like this, except with variable sizes.

Any help would be greatly appreciated! I have a lot to learn about matplotlib/pyplot and so on.

I have made a little bit of progress:

import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, PathPatch
from mpl_toolkits.mplot3d import Axes3D 
import mpl_toolkits.mplot3d.art3d as art3d

def cube(a,b,c,l): 
for zz in [c,c+l]: 
    for i in ["x","y","z"]:
        side = Rectangle((a, b), l,l)
        ax.add_patch(side)
        art3d.pathpatch_2d_to_3d(side, z=zz, zdir=i)


fig = plt.figure()
ax=fig.gca(projection='3d')

cube(0,0,0,1)



ax.set_xlim3d(-2, 2)
ax.set_ylim3d(-2, 2)
ax.set_zlim3d(-2, 2)

plt.show()

This plots a single cube.

EDIT:

More progress, I am now very close

import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, PathPatch
from mpl_toolkits.mplot3d import Axes3D 
import mpl_toolkits.mplot3d.art3d as art3d

cmap = plt.get_cmap('spring') #define the colors of the plot 
colors = [cmap(i) for i in np.linspace(0.1, 0.9, n+1)]  

def cube(a,b,c,l): #plots a cube of side l at (a,b,c)  
        for ll in [0,l]:
            for i in range(3):
                dire= ["x","y","z"]
                xdire = [b,a,a] 
                ydire = [c,c,b]
                zdire = [a,b,c]
                side = Rectangle((xdire[i],                                                         ydire[i]),facecolors[np.where(sizes == l)[0]],edgecolor='black')
                ax.add_patch(side)
                art3d.pathpatch_2d_to_3d(side, z=zdire[i]+ll, zdir=dire[i])

def plotter3D(X,Y,Z,sizes): #run cube(a,b,c,l) over the whole data set 
    for iX in range(len(X)):
        x = X[iX]
        y = Y[iX]
        z = Z[iX]
        for ix in range(len(x)): 
            cube(x[ix],y[ix],z[ix],sizes[iX])


fig = plt.figure() #open a figure 
ax=fig.gca(projection='3d') #make it 3d
plotter3D(X,Y,Z,sizes) #generate the cubes from the data set 
ax.set_xlim3d(0, length) #set the plot ranges 
ax.set_ylim3d(0, width)
ax.set_zlim3d(0, height)
plt.show()

This generates the desired output, although it seems to be see-through in some places when viewed from certain angles. You can see this in the small cube-- dead center at coordinates (1.5,2,3) Any idea how to fix this?

Another edit: The solution outined above has two problems: (1) I can't get equal aspect ratios for the three axes, and (2) The cubes are see-through from certain angles. Here's what the output looks like for a larger system