import itertools

vs = list(itertools.product(*[[-1,1]]*4))
es = [(v1, v2) for v1,v2 in itertools.product(vs,vs)
     if sum(1 for j in range(4) if v1[j]==v2[j])==3]


#A = (x0, y0, z0, 0)
#B = (x1, y1, 0, 0)
#C = (x2, 0, 0, 0)

x0,y0,z0,x1,y1,x2 = var('x0,y0,z0,x1,y1,x2')

f = (8 * ((1-x0)**2+(1-y0)**2+(1-z0)**2+1)**0.5
     +4* ((x1-x0)**2+(y1-y0)**2+z0**2)**0.5
     +2* ((x2-x1)**2+y1**2)**0.5
     + x2
    )

startP = tuple(random() for _ in range(6))
sol = minimize_constrained(f, [(0, 1)]*6, startP)
print (sol)
print (f(*sol))
#(0.8556652104160296, 0.7113316263361371, 0.42264981059562995, 0.7500064435495365, 0.5000068949778235, 0.49999972065448156)
#13.124355652992701





def getAllLines(paramVals):
    x00,y00,z00,x11,y11,x22 = paramVals
    yield ((x22,0,0,0), (-x22,0,0,0))
    for sgn3 in [1, -1]:
        for sgn2 in [1, -1]:
            yield ( (sgn3*x11, sgn2*y11,0,0), (sgn3*x22,0,0,0) )
            for sgn1 in [1, -1]:
                yield ((sgn1*x00, sgn2*y00, sgn3*z00,0),(sgn1*x11, sgn2*y11,0,0) )
                for sgn0 in [1, -1]:
                    yield ((sgn0*1,sgn1*1,sgn2*1,sgn3*1), (sgn0*x00, sgn1*y00, sgn2*z00,0) )

def project(p):
    u = 1.6
    v = 2.6
    t = 1
    cosu = float(cos(u))
    cosv = float(cos(v))
    cost = float(cos(t))
    sinu = float(sin(u))
    sinv = float(sin(v))
    sint = float(sin(t))
    return (
        cosu*cost*p[0] + cosu*sint*p[1] + sinu*cosv*p[2] + sinu*sinv*p[3],
        -sint*p[0] + cost*p[1],
        -sinu*cost*p[0] - sinu*sint*p[1] + cosu*cosv*p[2] + cosu*sinv*p[3],
        #-sinv*p[2] + cosv*p[3],
    )

#hcProj = polytopes.hypercube(4).schlegel_projection
#hcProj([1,1,1,-1]).plot().show()

g = Graphics()
g += points([project(v) for v in vs], size=20, color='blue')
for l in es: g += line([project(list(p)) for p in l], radius=0.005, color='blue' )
for l in getAllLines(sol): g += line([project(list(p)) for p in l], radius=0.02, color='green' )
g.show(frame=None)