1. import itertools
2.  
3. vs = list(itertools.product(*[[-1,1]]*4))
4. es = [(v1, v2) for v1,v2 in itertools.product(vs,vs)
5.      if sum(1 for j in range(4) if v1[j]==v2[j])==3]
6.  
7.  
8. #A = (x0, y0, z0, 0)
9. #B = (x1, y1, 0, 0)
10. #C = (x2, 0, 0, 0)
11.  
12. x0,y0,z0,x1,y1,x2 = var('x0,y0,z0,x1,y1,x2')
13.  
14. f = (8 * ((1-x0)**2+(1-y0)**2+(1-z0)**2+1)**0.5
15.      +4* ((x1-x0)**2+(y1-y0)**2+z0**2)**0.5
16.      +2* ((x2-x1)**2+y1**2)**0.5
17.      + x2
18.     )
19.  
20. startP = tuple(random() for _ in range(6))
21. sol = minimize_constrained(f, [(0, 1)]*6, startP)
22. print (sol)
23. print (f(*sol))
24. #(0.8556652104160296, 0.7113316263361371, 0.42264981059562995, 0.7500064435495365, 0.5000068949778235, 0.49999972065448156)
25. #13.124355652992701
26.  
27.  
28.  
29.  
30.  
31. def getAllLines(paramVals):
32.     x00,y00,z00,x11,y11,x22 = paramVals
33.     yield ((x22,0,0,0), (-x22,0,0,0))
34.     for sgn3 in [1, -1]:
35.         for sgn2 in [1, -1]:
36.             yield ( (sgn3*x11, sgn2*y11,0,0), (sgn3*x22,0,0,0) )
37.             for sgn1 in [1, -1]:
38.                 yield ((sgn1*x00, sgn2*y00, sgn3*z00,0),(sgn1*x11, sgn2*y11,0,0) )
39.                 for sgn0 in [1, -1]:
40.                     yield ((sgn0*1,sgn1*1,sgn2*1,sgn3*1), (sgn0*x00, sgn1*y00, sgn2*z00,0) )
41.  
42. def project(p):
43.     u = 1.6
44.     v = 2.6
45.     t = 1
46.     cosu = float(cos(u))
47.     cosv = float(cos(v))
48.     cost = float(cos(t))
49.     sinu = float(sin(u))
50.     sinv = float(sin(v))
51.     sint = float(sin(t))
52.     return (
53.         cosu*cost*p[0] + cosu*sint*p[1] + sinu*cosv*p[2] + sinu*sinv*p[3],
54.         -sint*p[0] + cost*p[1],
55.         -sinu*cost*p[0] - sinu*sint*p[1] + cosu*cosv*p[2] + cosu*sinv*p[3],
56.         #-sinv*p[2] + cosv*p[3],
57.     )
58.  
59. #hcProj = polytopes.hypercube(4).schlegel_projection
60. #hcProj([1,1,1,-1]).plot().show()
61.  
62. g = Graphics()
63. g += points([project(v) for v in vs], size=20, color='blue')
64. for l in es: g += line([project(list(p)) for p in l], radius=0.005, color='blue' )
65. for l in getAllLines(sol): g += line([project(list(p)) for p in l], radius=0.02, color='green' )
66. g.show(frame=None)
67.  
68.  
69.