
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import json


def draw_topo(topology,postTopo,labels,ax,node_color="#f2f0f0"):
    nx.draw(topology,postTopo,node_color=node_color,edge_color="#e3e3e3",node_size=400,ax=ax)
    nx.draw_networkx_labels(topology,postTopo,labels,font_color="#a0a0a0",ax=ax)
def draw_app(appG,pos,labels,ax):
    nx.draw(appG,pos,node_color="#4d997f", node_size=200,ax=ax)
    nx.draw_networkx_nodes(appG,pos,node_color="#4d997f",node_size=200,ax=ax)
    nx.draw_networkx_labels(appG,pos,labels,ax=ax)


# Infraestructure = PODs?
G = nx.balanced_tree(2,2)
labels = {l:l for l in G.nodes}
pos = nx.nx_agraph.graphviz_layout(G, prog="dot")
print(pos)
fig, ax = plt.subplots(figsize=(20,10))
draw_topo(G,pos,labels,ax)
plt.draw()

#Application: three services = GATEs
appG = nx.Graph()
nodes = {0:"frontend",1:"checkout",2:"ad"}
edges = [(0,1),(1,2)]
appG.add_nodes_from(nodes)
appG.add_edges_from(edges)

{0: (135.0, 162.0), 1: (99.0, 90.0), 2: (171.0, 90.0), 3: (27.0, 18.0), 4: (99.0, 18.0), 5: (171.0, 18.0), 6: (243.0, 18.0)}


# Quadratic model
# Ref: https://youtu.be/8O3So4Y2t_M?t=1456

# Connectivity Matrix
# Service 0 = with service 1
# service 1 = with service 2
# No weights
c = np.array([[0,1,0],[1,0,1],[0,1,0]])
c = nx.to_numpy_array(appG)

print(c)
print(c[:,0].sum())

[[0. 1. 0.]
 [1. 0. 1.]
 [0. 1. 0.]]
1.0


# Matriz
a = c*-1
print(a)

[[-0. -1. -0.]
 [-1. -0. -1.]
 [-0. -1. -0.]]


#Weights Gates to PADS aka. services to nodes 
# node3:service0 = 10
# node1:service1 = 5
# node4:service4 = 1
wires = [3,1,2]
# Los servicios han de estar conectados a un nodo? 
# opciones: 
    # comprobar todas las posibles opciones s0 node3, s0, node2, s0-nX 
    # ??

w = np.array([3,1,1])
diagonal = np.zeros(len(appG.nodes()))
for i in range(len(appG.nodes())):
    diagonal[i] = c[:,i].sum()+w[i]

np.fill_diagonal(a, diagonal)
print(a)

[[ 4. -1. -0.]
 [-1.  3. -1.]
 [-0. -1.  2.]]


# Bx vector (position)
# bx = w * xi
bx = np.zeros(len(appG.nodes()))
by = np.zeros(len(appG.nodes()))
for i in range(len(appG.nodes())):
    bx[i]=w[i]*pos[wires[i]][0]
    by[i]=w[i]*pos[wires[i]][1]

print(bx)
print(by)

[ 81.  99. 171.]
[54. 90. 90.]


x = np.linalg.solve(a, bx)
y = np.linalg.solve(a, by)
print(x)
print(y)
posApp = {i:(x[i],y[i]) for i in range(len(appG.nodes())) }
print(posApp)

[ 43.  91. 131.]
[30. 66. 78.]
{0: (43.0, 30.0), 1: (91.00000000000001, 66.0), 2: (131.00000000000003, 78.00000000000001)}


fig, ax = plt.subplots(figsize=(20,10))
draw_topo(G,pos,labels,ax)
draw_app(appG,posApp,ax)
plt.draw()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb Cell 10 in <cell line: 3>()
      <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#X11sZmlsZQ%3D%3D?line=0'>1</a> fig, ax = plt.subplots(figsize=(20,10))
      <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#X11sZmlsZQ%3D%3D?line=1'>2</a> draw_topo(G,pos,labels,ax)
----> <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#X11sZmlsZQ%3D%3D?line=2'>3</a> draw_app(appG,posApp,ax)
      <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#X11sZmlsZQ%3D%3D?line=3'>4</a> plt.draw()

TypeError: draw_app() missing 1 required positional argument: 'ax'


# MATRIX de WIRES
# Que ocurre si hay más conexiones entre alguno servicio y nodo?


# Weights Gates to PADS aka. services to nodes 
# node3:service0 = 10
# node1:service1 = 5
# node4:service4 = 1
wires = np.zeros(shape=(len(appG.nodes()),len(G.nodes())))
wires[0,3]=1
wires[0,2]=1
wires[0,1]=1
wires[1,5]=1
wires[2,0]=1
print(wires)

[[0. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0.]
 [1. 0. 0. 0. 0. 0. 0.]]


w = np.zeros(shape=(len(appG.nodes()),len(G.nodes())))
w[0,3]=100
w[0,2]=80
w[0,1]=180
w[1,5]=1
w[2,0]=1
print(w)

[[  0. 180.  80. 100.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   1.   0.]
 [  1.   0.   0.   0.   0.   0.   0.]]


c = np.array([[0,1,0],[1,0,1],[0,1,0]])
c = nx.to_numpy_array(appG)
a = c*-1
diagonal = np.zeros(len(appG.nodes()))
for i in range(len(appG.nodes())):
    diagonal[i] = c[:,i].sum()+w[i,:].sum()
np.fill_diagonal(a, diagonal)
print(diagonal)

bx = np.zeros(len(appG.nodes()))
by = np.zeros(len(appG.nodes()))

for i in range(len(appG.nodes())):
    print(i)
    sx,sy = 0,0
    for n in range(len(G.nodes())):
        sx += wires[i,n]*w[i,n]*pos[n][0]
        sy += wires[i,n]*w[i,n]*pos[n][1]
    
    bx[i]=sx
    by[i]=sy

[361.   3.   2.]
0
1
2


x = np.linalg.solve(a, bx)
y = np.linalg.solve(a, by)
print(x)
print(y)
posApp = {i:(x[i],y[i]) for i in range(len(appG.nodes())) }
print(posApp)

[ 95.10648918 133.44259567 134.22129784]
[ 69.99334443  67.59733777 114.79866889]
{0: (95.10648918469218, 69.99334442595674), 1: (133.4425956738769, 67.59733777038271), 2: (134.22129783693845, 114.79866888519135)}


wires[0,:].any()

True


def draw_forces(wires,pos,posApp,w,ax):
    apps, nodes = wires.shape
    for app in range(apps):
        for n in range(nodes):
            if wires[app,n]:
                x = [posApp[app][0],pos[n][0]]
                y = [posApp[app][1],pos[n][1]]
                ax.plot(x ,y, color="pink",linewidth=0.8,linestyle="--")
                ax.text(((x[1]+x[0])/2),((y[1]+y[0])/2),w[app,n],color="pink")


fig, ax = plt.subplots(figsize=(20,10))
draw_topo(G,pos,labels,ax)
draw_forces(wires,pos,posApp,w,ax)
draw_app(appG,posApp,ax)

plt.draw()


# Infraestructure = PODs?
G = nx.balanced_tree(2,2)
labels = {l:l for l in G.nodes}
pos = nx.circular_layout(G)
print(pos)
fig, ax = plt.subplots(figsize=(20,10))
draw_topo(G,pos,labels,ax)
plt.draw()

#Application: three services = GATEs
appG = nx.Graph()
nodes = {0:"frontend",1:"checkout",2:"ad"}
edges = [(0,1),(1,2)]
appG.add_nodes_from(nodes)
appG.add_edges_from(edges)

{0: array([1.00000000e+00, 8.51494896e-09]), 1: array([0.62348981, 0.78183148]), 2: array([-0.22252091,  0.97492788]), 3: array([-0.90096878,  0.43388381]), 4: array([-0.90096878, -0.43388373]), 5: array([-0.22252097, -0.97492786]), 6: array([ 0.62348963, -0.78183159])}


#Application: three services = GATEs
GI = nx.Graph()
nodes = {0:"A",1:"E",2:"G",3:"F",4:"C"}
lat = {(0,1):10,(1,2):5,(2,3):5,(3,4):10}
sumlat = sum(lat.values())
print(lat)

def f_acc_lat_IW(node1,node2):
    acclat = 0
    for x in range(node1,node2):
        acclat+=lat[(x,x+1)]
    return 40-acclat    

edges = [(node1,node2,{'weight': f_acc_lat_IW(node1,node2)}) for node1 in nodes for node2 in nodes if node1<node2]
print(edges)

# edges = [(0,1,{'weight': 10}),(1,2,{'weight': 5}),(2,3,{'weight': 5}),(3,4,{'weight': 10})]
# edges = [(0,1,{'weight': 10}),(1,2,{'weight': 5}),(2,3,{'weight': 5}),(3,4,{'weight': 10}),
        #  (0,2,{'weight': 15}),(0,3,{'weight': 20}),(0,4,{'weight': 30})
# ]

# fIW = lambda x:(x[0],x[1],{"weigth":40-x[2]["weight"]})
# edgesIW = list(map(fIW,edges))
# print(edgesIW)

{(0, 1): 10, (1, 2): 5, (2, 3): 5, (3, 4): 10}
[(0, 1, {'weight': 30}), (0, 2, {'weight': 25}), (0, 3, {'weight': 20}), (0, 4, {'weight': 10}), (1, 2, {'weight': 35}), (1, 3, {'weight': 30}), (1, 4, {'weight': 20}), (2, 3, {'weight': 35}), (2, 4, {'weight': 25}), (3, 4, {'weight': 30})]


GI.add_nodes_from(nodes)
GI.add_edges_from(edges)
posR = nx.circular_layout(GI)
pos = nx.spring_layout(GI,pos=posR,iterations=100) 

fig, ax = plt.subplots(figsize=(20,10))
draw_topo(GI,pos,nodes,ax)
plt.draw()


# MATRIX de WIRES
# Que ocurre si hay más conexiones entre alguno servicio y nodo?
#Application: three services = GATEs
appG = nx.Graph()
nodesG = {0:"frontend",1:"backend"}
edgesG = [(0,1)]
appG.add_nodes_from(nodesG)
appG.add_edges_from(edgesG)


# Weights Gates to PADS aka. services to nodes 
# node3:service0 = 10
# node1:service1 = 5
# node4:service4 = 1
wires = np.zeros(shape=(len(appG.nodes()),len(GI.nodes())))
wires[0,0]=1 #at A
wires[0,4]=1 #at C
wires[1,2]=1 #at G
print(wires)

[[1. 0. 0. 0. 1.]
 [0. 0. 1. 0. 0.]]


w = np.zeros(shape=(len(appG.nodes()),len(GI.nodes())))
w[0,0]=100
w[0,4]=100
w[1,3]=100
print(w)

[[100.   0.   0.   0. 100.]
 [  0.   0.   0. 100.   0.]]


c = np.array([[0,1],[1,0]]) # matriz de connexiones

c = nx.to_numpy_array(appG)
a = c*-1
diagonal = np.zeros(len(appG.nodes()))
for i in range(len(appG.nodes())):
    diagonal[i] = c[:,i].sum()+w[i,:].sum()
np.fill_diagonal(a, diagonal)
print(diagonal)

bx = np.zeros(len(appG.nodes()))
by = np.zeros(len(appG.nodes()))

for i in range(len(appG.nodes())):
    print(i)
    sx,sy = 0,0
    for n in range(len(GI.nodes())):
        sx += wires[i,n]*w[i,n]*pos[n][0]
        sy += wires[i,n]*w[i,n]*pos[n][1]
    
    bx[i]=sx
    by[i]=sy

[201. 101.]
0
1


x = np.linalg.solve(a, bx)
y = np.linalg.solve(a, by)
print(x)
print(y)
posApp = {i:(x[i],y[i]) for i in range(len(appG.nodes())) }
print(posApp)

[-0.20726364 -0.00205212]
[0.15053179 0.00149041]
{0: (-0.20726363871567047, 0.1505317882330052), 1: (-0.0020521152348086185, 0.0014904137448812393)}


fig, ax = plt.subplots(figsize=(20,10))
draw_topo(GI,pos,nodes,ax)
draw_forces(wires,pos,posApp,w,ax)
draw_app(appG,posApp,nodesG,ax)
plt.draw()


#Application: three services = GATEs
GI = nx.Graph()
nodes = {0:"E",1:"G",2:"F"}
lat = {(0,1):5,(1,2):5}
sumlat = sum(lat.values())

def f_acc_lat_IW(node1,node2):
    acclat = 0
    for x in range(node1,node2):
        acclat+=lat[(x,x+1)]
    return sumlat-acclat    

edges = [(node1,node2,{'weight': f_acc_lat_IW(node1,node2)}) for node1 in nodes for node2 in nodes if node1<node2]
print(edges)
GI.add_nodes_from(nodes)
GI.add_edges_from(edges)
posR = nx.circular_layout(GI)
pos = nx.spring_layout(GI,pos=posR,iterations=100) 

fig, ax = plt.subplots(figsize=(20,10))
draw_topo(GI,pos,nodes,ax)
plt.draw()

appG = nx.Graph()
nodesG = {0:"frontend",1:"backend",2:"Datastore"}
edgesG = [(0,1),(1,2)]
appG.add_nodes_from(nodesG)
appG.add_edges_from(edgesG)

wires = np.zeros(shape=(len(appG.nodes()),len(GI.nodes())))
# wires[0,0]=1 #at E
# wires[0,2]=1 #at F
wires[0,1]=1 #at F
wires[1,1]=1 #at G
wires[2,2]=1 #at G
w = np.zeros(shape=(len(appG.nodes()),len(GI.nodes())))
# w[0,0]=80
# w[0,2]=60
w[0,1]=20
w[1,1]=100
w[2,2]=200
c = np.array([[0,1,0],[1,0,1],[0,0,1]]) # matriz de connexiones

c = nx.to_numpy_array(appG)
a = c*-1
diagonal = np.zeros(len(appG.nodes()))
for i in range(len(appG.nodes())):
    diagonal[i] = c[:,i].sum()+w[i,:].sum()
np.fill_diagonal(a, diagonal)
print(diagonal)

bx = np.zeros(len(appG.nodes()))
by = np.zeros(len(appG.nodes()))

for i in range(len(appG.nodes())):
    print(i)
    sx,sy = 0,0
    for n in range(len(GI.nodes())):
        sx += wires[i,n]*w[i,n]*pos[n][0]
        sy += wires[i,n]*w[i,n]*pos[n][1]
    
    bx[i]=sx
    by[i]=sy

x = np.linalg.solve(a, bx)
y = np.linalg.solve(a, by)
posApp = {i:(x[i],y[i]) for i in range(len(appG.nodes())) }

fig, ax = plt.subplots(figsize=(20,10))
draw_topo(GI,pos,nodes,ax)
draw_forces(wires,pos,posApp,w,ax)
draw_app(appG,posApp,nodesG,ax)
plt.draw()

[(0, 1, {'weight': 5}), (0, 2, {'weight': 0}), (1, 2, {'weight': 5})]
[ 21. 102. 201.]
0
1
2


import sympy as sp


# Infraestructure = PODs?
G = nx.balanced_tree(2,2)
labels = {l:l for l in G.nodes}
pos = nx.nx_agraph.graphviz_layout(G, prog="dot")
print(pos)
fig, ax = plt.subplots(figsize=(20,10))
draw_topo(G,pos,labels,ax)
plt.draw()

{0: (135.0, 162.0), 1: (99.0, 90.0), 2: (171.0, 90.0), 3: (27.0, 18.0), 4: (99.0, 18.0), 5: (171.0, 18.0), 6: (243.0, 18.0)}


# latency matrix
paths = dict(nx.all_pairs_shortest_path(G))

lats = [len(paths[n1][n2])-1 for n1 in G.nodes() for n2 in G.nodes()]
latM = np.array(lats).reshape(len(G.nodes()),len(G.nodes()))
latM = np.triu(latM)

print(latM)

[[0 1 1 2 2 2 2]
 [0 0 2 1 1 3 3]
 [0 0 0 3 3 1 1]
 [0 0 0 0 2 4 4]
 [0 0 0 0 0 4 4]
 [0 0 0 0 0 0 2]
 [0 0 0 0 0 0 0]]


# Solve euclidean equations 
dim = 2

vars = "abcdefghijklmnopqrstuvwxyz"
if len(G.nodes())>len(vars):
    assert("Warning. add more names ")

unkX = [list(vars)[e][0]+"%i"%d for e,_ in enumerate(G.nodes()) for d in range(dim) ]
unkX

['a0',
 'a1',
 'b0',
 'b1',
 'c0',
 'c1',
 'd0',
 'd1',
 'e0',
 'e1',
 'f0',
 'f1',
 'g0',
 'g1']


syms = sp.symbols(unkX)
equations = []
for ix, row in enumerate(latM):
    for iy, distance in enumerate(row):
        if iy <= ix: continue
        if distance > 0:
            sumDim = 0
            for d in range(dim):
                sumDim += (syms[(ix*dim)+d] - syms[(iy*dim)+d])**2
            eq = sp.sqrt(sumDim)
            print(eq,distance)
            equations.append(sp.Eq(eq, distance))

        ##### WRONG EQUATION
        # if distance > 0:
        #     for d in range(dim):
        #         eq = syms[(ix*dim)+d]**2 - 2 * syms[(ix*dim)+d] * syms[(iy*dim)+d] - syms[(iy*dim)+d]**2
        #         print(eq,distance)
        #         equations.append(sp.Eq(eq, distance))

sqrt((a0 - b0)**2 + (a1 - b1)**2) 1
sqrt((a0 - c0)**2 + (a1 - c1)**2) 1
sqrt((a0 - d0)**2 + (a1 - d1)**2) 2
sqrt((a0 - e0)**2 + (a1 - e1)**2) 2
sqrt((a0 - f0)**2 + (a1 - f1)**2) 2
sqrt((a0 - g0)**2 + (a1 - g1)**2) 2
sqrt((b0 - c0)**2 + (b1 - c1)**2) 2
sqrt((b0 - d0)**2 + (b1 - d1)**2) 1
sqrt((b0 - e0)**2 + (b1 - e1)**2) 1
sqrt((b0 - f0)**2 + (b1 - f1)**2) 3
sqrt((b0 - g0)**2 + (b1 - g1)**2) 3
sqrt((c0 - d0)**2 + (c1 - d1)**2) 3
sqrt((c0 - e0)**2 + (c1 - e1)**2) 3
sqrt((c0 - f0)**2 + (c1 - f1)**2) 1
sqrt((c0 - g0)**2 + (c1 - g1)**2) 1
sqrt((d0 - e0)**2 + (d1 - e1)**2) 2
sqrt((d0 - f0)**2 + (d1 - f1)**2) 4
sqrt((d0 - g0)**2 + (d1 - g1)**2) 4
sqrt((e0 - f0)**2 + (e1 - f1)**2) 4
sqrt((e0 - g0)**2 + (e1 - g1)**2) 4
sqrt((f0 - g0)**2 + (f1 - g1)**2) 2


ans = sp.solve(equations, syms)
ans

[]


#Application: three services = GATEs
GI = nx.Graph()
nodes = {0:"A",1:"B",2:"C"}
edges = [(0,1,{'weight': 1}),(1,2,{'weight': 1}),(2,0,{'weight': 1})]
GI.add_nodes_from(nodes)
GI.add_edges_from(edges)
labels = {l:l for l in GI.nodes}
pos = nx.nx_agraph.graphviz_layout(GI, prog="dot")

fig, ax = plt.subplots(figsize=(10,6))
draw_topo(GI,pos,labels,ax)
plt.draw()


# latency matrix
paths = dict(nx.all_pairs_shortest_path(GI))

lats = [len(paths[n1][n2])-1 for n1 in GI.nodes() for n2 in GI.nodes()]
latM = np.array(lats).reshape(len(GI.nodes()),len(GI.nodes()))
latM = np.triu(latM)

print(latM)

# Solve euclidean equations 
dim = 2

vars = "abcdefghijklmnopqrstuvwxyz"
if len(GI.nodes())>len(vars):
    assert("Warning. add more names ")

unkX = [list(vars)[e][0]+"%i"%d for e,_ in enumerate(GI.nodes()) for d in range(dim) ]
print(unkX)

[[0 1 1]
 [0 0 1]
 [0 0 0]]
['a0', 'a1', 'b0', 'b1', 'c0', 'c1']


syms = sp.symbols(unkX,real=True)
equations = []
for ix, row in enumerate(latM):
    for iy, distance in enumerate(row):
        if iy <= ix: continue
        if distance > 0:
            sumDim = 0
            for d in range(dim):
                sumDim += (syms[(ix*dim)+d] - syms[(iy*dim)+d])**2
            eq = sp.sqrt(sumDim)
            print(eq,distance)
            equations.append(sp.Eq(eq, distance))

sqrt((a0 - b0)**2 + (a1 - b1)**2) 1
sqrt((a0 - c0)**2 + (a1 - c1)**2) 1
sqrt((b0 - c0)**2 + (b1 - c1)**2) 1


ans = sp.solve(equations, syms)
print("Number of solutions ",len(ans))

Number of solutions  4


ndim = 2

for sol in range(len(ans)):
    assignments = {}
    values = 0
    for ix,solsymbol in enumerate(ans[sol]):
        if type(solsymbol)== sp.Symbol:
            assignments[syms[ix]]=values 
            if ix%ndim != 0:
                values+=1
    print(assignments)

    points = tuple(float(t.subs(assignments).evalf(4)) for t in ans[sol])
    print(points)

    pos2 = {}
    for ix,n in enumerate(GI.nodes()):
        t = ()
        for d in range(ndim):
            t = t+(points[(ix*ndim)+d],)  
        pos2[n]=t
            
    print(pos2)
    break

{b1: 0, c0: 1, c1: 1}
(1.86602783203125, 0.5, 1.0, 0.0, 1.0, 1.0)
{0: (1.86602783203125, 0.5), 1: (1.0, 0.0), 2: (1.0, 1.0)}


fig, ax = plt.subplots(figsize=(10,6))
draw_topo(GI,pos2,labels,ax,"#ff00f0")
ax.set_aspect('equal')
ax.grid(True, which='both')
ax.axhline(y=0, color='k')
ax.axvline(x=1, color='k')
ax.tick_params(left=True, bottom=True, labelleft=True, labelbottom=True)
plt.axis('on')
plt.draw()


import time


# Infraestructure = PODs?
G = nx.balanced_tree(2,2)
labels = {l:l for l in G.nodes}
# latency matrix
paths = dict(nx.all_pairs_shortest_path(G))
lats = [len(paths[n1][n2])-1 for n1 in G.nodes() for n2 in G.nodes()]
latM = np.array(lats).reshape(len(G.nodes()),len(G.nodes()))
latM = np.triu(latM)

ansG = None
# Solve euclidean equations 
for dim in range(2,10):
# dim = 2
    print("Dimension: %i"%dim)
    vars = "abcdefghijklmnopqrstuvwxyz"
    if len(G.nodes())>len(vars):
        assert("Warning. add more names ")

    unkX = [list(vars)[e][0]+"%i"%d for e,_ in enumerate(G.nodes()) for d in range(dim) ]
    syms = sp.symbols(unkX,real=True)
    equations = []
    for ix, row in enumerate(latM):
        for iy, distance in enumerate(row):
            if iy <= ix: continue
            if distance > 0:
                sumDim = 0
                for d in range(dim):
                    sumDim += (syms[(ix*dim)+d] - syms[(iy*dim)+d])**2
                eq = sp.sqrt(sumDim)
                equations.append(sp.Eq(eq, distance))

    start = time.time()
    ans = sp.solve(equations, syms,particular = True, quick=True,warn=True  )
    print("\tTime: ",time.time()-start)
    if len(ans) == 0:
        print("\tNo solution for a dim == %i"%dim)
    else:
        ansG = ans
        break

Dimension: 2
	Time:  0.15458893775939941
	No solution for a dim == 2
Dimension: 3
	Time:  26.652655839920044
	No solution for a dim == 3
Dimension: 4

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb Cell 46 in <cell line: 12>()
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y105sZmlsZQ%3D%3D?line=29'>30</a>             equations.append(sp.Eq(eq, distance))
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y105sZmlsZQ%3D%3D?line=31'>32</a> start = time.time()
---> <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y105sZmlsZQ%3D%3D?line=32'>33</a> ans = sp.solve(equations, syms,particular = True, quick=True,warn=True  )
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y105sZmlsZQ%3D%3D?line=33'>34</a> print("\tTime: ",time.time()-start)
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y105sZmlsZQ%3D%3D?line=34'>35</a> if len(ans) == 0:

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/solvers/solvers.py:1114, in solve(f, *symbols, **flags)
   1112     solution = _solve(f[0], *symbols, **flags)
   1113 else:
-> 1114     solution = _solve_system(f, symbols, **flags)
   1116 #
   1117 # postprocessing
   1118 ###########################################################################
   1119 # Restore masked-off objects
   1120 if non_inverts:

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/solvers/solvers.py:1868, in _solve_system(exprs, symbols, **flags)
   1865 for syms in subsets(free, len(polys)):
   1866     try:
   1867         # returns [] or list of tuples of solutions for syms
-> 1868         res = solve_poly_system(polys, *syms)
   1869         if res:
   1870             for r in res:

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/solvers/polysys.py:76, in solve_poly_system(seq, strict, *gens, **args)
     73         except SolveFailed:
     74             pass
---> 76 return solve_generic(polys, opt, strict=strict)

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/solvers/polysys.py:326, in solve_generic(polys, opt, strict)
    323     return solutions
    325 try:
--> 326     result = _solve_reduced_system(polys, opt.gens, entry=True)
    327 except CoercionFailed:
    328     raise NotImplementedError

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/solvers/polysys.py:266, in solve_generic.<locals>._solve_reduced_system(system, gens, entry)
    263     zeros = list(roots(system[0], gens[-1], strict=strict).keys())
    264     return [(zero,) for zero in zeros]
--> 266 basis = groebner(system, gens, polys=True)
    268 if len(basis) == 1 and basis[0].is_ground:
    269     if not entry:

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/polytools.py:6922, in groebner(F, *gens, **args)
   6871 @public
   6872 def groebner(F, *gens, **args):
   6873     """
   6874     Computes the reduced Groebner basis for a set of polynomials.
   6875 
   (...)
   6920 
   6921     """
-> 6922     return GroebnerBasis(F, *gens, **args)

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/polytools.py:6961, in GroebnerBasis.__new__(cls, F, *gens, **args)
   6957 ring = PolyRing(opt.gens, opt.domain, opt.order)
   6959 polys = [ring.from_dict(poly.rep.to_dict()) for poly in polys if poly]
-> 6961 G = _groebner(polys, ring, method=opt.method)
   6962 G = [Poly._from_dict(g, opt) for g in G]
   6964 return cls._new(G, opt)

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/groebnertools.py:43, in groebner(seq, ring, method)
     40     else:
     41         seq = [ s.set_ring(ring) for s in seq ]
---> 43 G = _groebner(seq, ring)
     45 if orig is not None:
     46     G = [ g.clear_denoms()[1].set_ring(orig) for g in G ]

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/groebnertools.py:238, in _buchberger(f, ring)
    236 # ordering divisors is on average more efficient [Cox] page 111
    237 G1 = sorted(G, key=lambda g: order(f[g].LM))
--> 238 ht = normal(h, G1)
    240 if ht:
    241     G, CP = update(G, CP, ht[1])

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/groebnertools.py:104, in _buchberger.<locals>.normal(g, J)
    103 def normal(g, J):
--> 104     h = g.rem([ f[j] for j in J ])
    106     if not h:
    107         return None

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/rings.py:1493, in PolyElement.rem(self, G)
   1491 for mg, cg in g.iterterms():
   1492     m1 = monomial_mul(mg, m)
-> 1493     c1 = get(m1, zero) - c*cg
   1494     if not c1:
   1495         del f[m1]

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/fields.py:451, in FracElement.__sub__(f, g)
    449         return f.new(f.numer - g.numer, f.denom)
    450     else:
--> 451         return f.new(f.numer*g.denom - f.denom*g.numer, f.denom*g.denom)
    452 elif isinstance(g, field.ring.dtype):
    453     return f.new(f.numer - f.denom*g, f.denom)

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/fields.py:301, in FracElement.new(f, numer, denom)
    300 def new(f, numer, denom):
--> 301     return f.raw_new(*numer.cancel(denom))

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/rings.py:2223, in PolyElement.cancel(self, g)
   2220 domain = ring.domain
   2222 if not (domain.is_Field and domain.has_assoc_Ring):
-> 2223     _, p, q = f.cofactors(g)
   2224 else:
   2225     new_ring = ring.clone(domain=domain.get_ring())

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/rings.py:2138, in PolyElement.cofactors(f, g)
   2135     h, cfg, cff = g._gcd_monom(f)
   2136     return h, cff, cfg
-> 2138 J, (f, g) = f.deflate(g)
   2139 h, cff, cfg = f._gcd(g)
   2141 return (h.inflate(J), cff.inflate(J), cfg.inflate(J))

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/sympy/polys/rings.py:2069, in PolyElement.deflate(f, *G)
   2067     for monom in p.itermonoms():
   2068         for i, m in enumerate(monom):
-> 2069             J[i] = igcd(J[i], m)
   2071 for i, b in enumerate(J):
   2072     if not b:

KeyboardInterrupt:


from scipy.optimize import fsolve
import time
def equations(vars):
    a0,a1,b0,b1,c0,c1 = vars
    equations = np.empty((3))
    equations[0] = np.sqrt((a0 - b0)**2 + (a1 - b1)**2)-1
    equations[1] = np.sqrt((a0 - c0)**2 + (a1 - c1)**2) - 1
    equations[2] = np.sqrt((b0 - c0)**2 + (b1 - c1)**2) -2
    return equations

ini = (1, 1, 1,1,1,1)
x, y, z =  fsolve(equations, ini)
print(x, y, z)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb Cell 48 in <cell line: 12>()
      <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y111sZmlsZQ%3D%3D?line=8'>9</a>     return equations
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y111sZmlsZQ%3D%3D?line=10'>11</a> ini = (1, 1, 1,1,1,1)
---> <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y111sZmlsZQ%3D%3D?line=11'>12</a> x, y, z =  fsolve(equations, ini)
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y111sZmlsZQ%3D%3D?line=12'>13</a> print(x, y, z)

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:160, in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
     49 """
     50 Find the roots of a function.
     51 
   (...)
    150 
    151 """
    152 options = {'col_deriv': col_deriv,
    153            'xtol': xtol,
    154            'maxfev': maxfev,
   (...)
    157            'factor': factor,
    158            'diag': diag}
--> 160 res = _root_hybr(func, x0, args, jac=fprime, **options)
    161 if full_output:
    162     x = res['x']

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:226, in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    224 if not isinstance(args, tuple):
    225     args = (args,)
--> 226 shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
    227 if epsfcn is None:
    228     epsfcn = finfo(dtype).eps

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:38, in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     36             msg += "."
     37         msg += 'Shape should be %s but it is %s.' % (output_shape, shape(res))
---> 38         raise TypeError(msg)
     39 if issubdtype(res.dtype, inexact):
     40     dt = res.dtype

TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'equations'.Shape should be (6,) but it is (3,).


from gekko import GEKKO

m = GEKKO()
x = m.Var(value=0)      # define new variable, initial value=0
y = m.Var(value=1)      # define new variable, initial value=1
z = m.Var(value=1)
m.Equations([x + 2*y==0, x**2+y**2==1, z==y**2]) # equations
m.solve(disp=False)     # solve

print([x.value[0],y.value[0]]) # print solution

[-0.89442724267, 0.44721362133]


# # Infraestructure = PODs?
# G = nx.balanced_tree(2,2)
# labels = {l:l for l in G.nodes}
# # latency matrix
# paths = dict(nx.all_pairs_shortest_path(G))
# lats = [len(paths[n1][n2])-1 for n1 in G.nodes() for n2 in G.nodes()]
# latM = np.array(lats).reshape(len(G.nodes()),len(G.nodes()))
# latM = np.triu(latM)


#Application: three services = GATEs
G = nx.Graph()
nodes = {0:"A",1:"B",2:"C"}
edges = [(0,1,{'weight': 1}),(1,2,{'weight': 1}),(2,0,{'weight': 1})]
G.add_nodes_from(nodes)
G.add_edges_from(edges)
labels = {l:l for l in G.nodes}
pos = nx.nx_agraph.graphviz_layout(G, prog="dot")


# latency matrix
paths = dict(nx.all_pairs_shortest_path(G))

lats = [len(paths[n1][n2])-1 for n1 in G.nodes() for n2 in G.nodes()]
latM = np.array(lats).reshape(len(G.nodes()),len(G.nodes()))
latM = np.triu(latM)


m = GEKKO()

# Solve euclidean equations 
for dim in range(2,10):
# dim = 2
    print("Dimension: %i"%dim)
    vars = "abcdefghijklmnopqrstuvwxyz"
    if len(G.nodes())>len(vars):
        assert("Warning. add more names ")

    unkX = [list(vars)[e][0]+"%i"%d for e,_ in enumerate(G.nodes()) for d in range(dim) ]
    syms = [m.Var(lb=-100, ub=100) for i in range(len(unkX))]

    print(syms)

    equations = []
    for ix, row in enumerate(latM):
        for iy, distance in enumerate(row):
            if iy <= ix: continue
            if distance > 0:
                sumDim = 0
                for d in range(dim):
                    sumDim += (syms[(ix*dim)+d] - syms[(iy*dim)+d]) #**2
                # eq = sp.sqrt(sumDim)
                eq = sumDim
                m.Equation(eq==distance)
                # equations.append(eq == distance)

    
    print(equations)
    m.solve(disp=False)
    for var in syms:
        print(var.value)

    break

m.cleanup()

Dimension: 2
[0, 0, 0, 0, 0, 0]
[]

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb Cell 51 in <cell line: 32>()
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y116sZmlsZQ%3D%3D?line=54'>55</a>             # equations.append(eq == distance)
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y116sZmlsZQ%3D%3D?line=57'>58</a> print(equations)
---> <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y116sZmlsZQ%3D%3D?line=58'>59</a> m.solve(disp=False)
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y116sZmlsZQ%3D%3D?line=59'>60</a> for var in syms:
     <a href='vscode-notebook-cell:/Users/isaac/Projects/ForceDirectAllocation/QuadraticPlacement.ipynb#Y116sZmlsZQ%3D%3D?line=60'>61</a>     print(var.value)

File ~/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages/gekko/gekko.py:2185, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
   2183 #print APM error message and die
   2184 if (debug >= 1) and ('@error' in response):
-> 2185     raise Exception(response)
   2187 #load results
   2188 def byte2str(byte):

Exception:  @error: Solution Not Found


m.cleanup()

<bound method GEKKO.Equations of <gekko.gekko.GEKKO object at 0x16b80f5e0>>

Quadratic Placement for service allocation problem¶

Multiples wires¶

Test: Circular layout¶

Equidistance of latencies as a function of topology¶

Solving undetermined number of dimensions¶

Testing previous code: equilateral triangle¶

Scaling the problem: Dimension up!¶

Other nonlinear solutions ?¶

GEKKO Suite¶