#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Matthew Donoughe

from FuncDesigner import *
from openopt import MILP
from myooutil import *
import copy
import numpy

# network size, in nodes
NET_SIZE = 6
# number of wavelengths per link
WAVELENGTHS = 1

# 0 -- 1 -- 2
# | \  |  / |
# |  \ | /  |
# 3 -- 4 -- 5

# create the physical network matrix
# multiple links from s to d are allowed
# asymetric networks are allowed
# the diagonal is ignored
physical = [
        [0, 1, 0, 1, 1, 0], # 0 to d
        [1, 0, 1, 0, 1, 0], # 1 to d
        [0, 1, 0, 0, 1, 1], # 2 to d
        [1, 0, 0, 0, 1, 0], # 3 to d
        [1, 1, 1, 1, 0, 1], # 4 to d
        [0, 0, 1, 0, 1, 0]  # 5 to d
        ]
print('physical connections:')
print_matrix(physical)

# create the demand network matrix
# multiple requests from s to d are allowed
# asymmetric networks are allowed
# the diagonal is ignored(nodes are assumed to have an internal self-connection)
demand = [
        [0, 0, 0, 0, 0, 1],
        [0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0],
        [1, 0, 0, 0, 0, 0],
        ]
# it is also possible to define the demand matrix like this:
#demand = make_network(NET_SIZE, (0, 5), (0, 5), (0, 5), (5, 0), (0, 3), (0, 3), (2, 4), (4, 2), (0, 4), (0, 4), (4, 0))
print('demand:')
print_matrix(demand)

constraints = []

# allocate space for all variables in a single vector
# flow_max is one element
# connection_wavelengths is NET_SIZE ** 2 * WAVELENGTHS logical elements
# NETSIZE places where WAVELENGTHS positions have the same s and d
# links is NET_SIZE ** 4 * WAVELENGTHS logical elements
# NET_SIZE places where NET_SIZE ** 2 * WAVELENGTHS positions have the same s and d
# NET_SIZE ** 2 * WAVELENGTHS places where NETSIZE positions have the same i and j
# NET_SIZE ** 2 * WAVELENGTHS places where 1 position has the same s d i j(and was subtracted twice)
oodata = oovar(size = NET_SIZE ** 4 * WAVELENGTHS - 2 * NET_SIZE ** 3 * WAVELENGTHS + 2 * NET_SIZE ** 2 * WAVELENGTHS - NET_SIZE * WAVELENGTHS + 1)
#oodata = [oovar('oodata[%d]' % i) for i in xrange(NET_SIZE ** 4 * WAVELENGTHS - 2 * NET_SIZE ** 3 * WAVELENGTHS + 2 * NET_SIZE ** 2 * WAVELENGTHS - NET_SIZE * WAVELENGTHS + 1)]
data_index = 0

# assign wavelengnths to connections
connection_wavelengths = make_matrix(lambda s, d, w: 0, NET_SIZE, NET_SIZE, WAVELENGTHS)
# strip out self connections
for s in range(NET_SIZE):
    connection_wavelengths[s][s] = []
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                connection_wavelengths[s][d][w] = oodata[data_index]
                data_index += 1

# assign connections to physical links
links = make_matrix(lambda s, d, w, i, j: 0, NET_SIZE, NET_SIZE, WAVELENGTHS, NET_SIZE, NET_SIZE)
# strip out self connections
for s in range(NET_SIZE):
    links[s][s] = make_matrix(lambda w, i: [], WAVELENGTHS, NET_SIZE)
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                for i in range(NET_SIZE):
                    links[s][d][w][i][i] = 0
                    for j in range(NET_SIZE):
                        if i != j:
                            links[s][d][w][i][j] = oodata[data_index]
                            data_index += 1

# 1
flow_max = oodata[data_index]
data_index += 1
print(data_index)

# 2
# Fmax ≥ total number of flows in the same direction on any link
temp = make_matrix(lambda i, j: [], NET_SIZE, NET_SIZE)
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                for i in range(NET_SIZE):
                    for j in range(NET_SIZE):
                        if i != j:
                            temp[i][j] += [links[s][d][w][i][j]]
for i in range(NET_SIZE):
    for j in range(NET_SIZE):
        if i != j:
            constraints += [(flow_max >= sum(temp[i][j]))('Fmax >= flows on %d->%d' % (i, j))]
del temp

# 3
# the difference between inbound and outbound must be the number of connections added/removed to/from the network at that node
inbound = make_matrix(lambda s, d, w, j: [], NET_SIZE, NET_SIZE, WAVELENGTHS, NET_SIZE)
outbound = make_matrix(lambda s, d, w, i: [], NET_SIZE, NET_SIZE, WAVELENGTHS, NET_SIZE)
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                for i in range(NET_SIZE):
                    for j in range(NET_SIZE):
                        if i != j:
                            inbound[s][d][w][j] += [links[s][d][w][i][j]]
                            outbound[s][d][w][i] += [links[s][d][w][i][j]]
inbound = [[[[sum(j) for j in w] for w in d] for d in s] for s in inbound]
outbound = [[[[sum(i) for i in w] for w in d] for d in s] for s in outbound]
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                for j in range(NET_SIZE):
                    diff = (inbound[s][d][w][j] - outbound[s][d][w][j])('diff[s=%d,d=%d,w=%d,j=%d]' % (s, d, w, j))
                    if s == j:
                        constraints += [(diff == -connection_wavelengths[s][d][w])('%s == -%s' % (diff, connection_wavelengths[s][d][w]))]
                    elif d == j:
                        constraints += [(diff == connection_wavelengths[s][d][w])('%s == %s' % (diff, connection_wavelengths[s][d][w]))]
                    else:
                        constraints += [(diff == 0)('%s == 0' % diff)]
del inbound
del outbound

# 4
# wavelengths used between s and d must equal the demand for connections
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            constraints += [(sum(connection_wavelengths[s][d]) == demand[s][d])('connections %d->%d == demand' % (s, d))]

# 5
# connections from i to j on wavelength w must not be negative
# in the homework formulation, this also says that the number must be <= 1,
# but a) this is already taken care of by ensuring only one lightpath per wavelength per link,
# and b) this is incorrect if there are multiple physical links(then there could be more)
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                for i in range(NET_SIZE):
                    for j in range(NET_SIZE):
                        if i != j:
                            constraints += [(links[s][d][w][i][j] >= 0)('%s >= 0' % links[s][d][w][i][j])]

# 6
# number of lightpaths on a given wavelength may not exceed the number of links
# this differs from the homework formulation in two ways:
# 1) it prevents lightpaths over links that do not exist
# 2) it allows for multiple physical links
temp = make_matrix(lambda i, j, w: [], NET_SIZE, NET_SIZE, WAVELENGTHS)
for s in range(NET_SIZE):
    for d in range(NET_SIZE):
        if s != d:
            for w in range(WAVELENGTHS):
                for i in range(NET_SIZE):
                    for j in range(NET_SIZE):
                        if i != j:
                            temp[i][j][w] += [links[s][d][w][i][j]]
for i in range(NET_SIZE):
    for j in range(NET_SIZE):
        if i != j:
            for w in range(WAVELENGTHS):
                constraints += [(sum(temp[i][j][w]) <= physical[i][j])('lightpaths[i=%d,j=%d,w=%d] <= physical[i=%d,j=%d]' % (i, j, w, i, j))]
del temp

# start with all variables set to 0
startPoint = {oodata: numpy.zeros(data_index)}
#startPoint = dict([(v, 0) for v in oodata])

print('OpenOpt is converting %d constraints and %d variables...' % (len(constraints), len(startPoint)))
problem = MILP(flow_max, startPoint, constraints = constraints, intVars = [oodata])
#problem = MILP(flow_max, startPoint, constraints = constraints, intVars = oodata)
# cplex requires:
# OpenOpt 0.33
# CPLEX 12
# the Python API distributed with CPLEX must be installed(on Mac OS X x86_64, you must run 'arch -i386 python2.6 solve.py' or it won't work)
# the CPLEX license must be installed properly
result = problem.minimize('cplex')
# glpk requires:
# OpenOpt
# CVXOPT
# glpk
#result = problem.minimize('glpk')

# extract the results from OpenOpt
links = clean_matrix(links, result)
# convert the results to integers
links = [[[[[int(j) for j in i] for i in w] for w in d] for d in s] for s in links]

# graph the result
print('writing graph to graph.dot...')
# make copies so we can better draw bidirectional links without side effects
phys = copy.deepcopy(physical)
lnks = copy.deepcopy(links)
with open('graph.dot', 'w') as g:
    g.write('digraph {\n')
    g.write('edge [colorscheme=set19,decorate=true,fontsize=6];\n')
    g.write('n0 [pos="0,3",pin=true];\n')
    g.write('n1 [pos="3,3",pin=true];\n')
    g.write('n2 [pos="6,3",pin=true];\n')
    g.write('n3 [pos="0,0",pin=true];\n')
    g.write('n4 [pos="3,0",pin=true];\n')
    g.write('n5 [pos="6,0",pin=true];\n')

    # draw physical links
    g.write('// physical links\n')
    for s in range(NET_SIZE):
        for d in range(NET_SIZE):
            if s != d:
                for i in range(phys[s][d]):
                    if phys[d][s] > 0:
                        g.write('n%d -> n%d [dir=none];\n' % (s, d))
                        phys[d][s] -= 1
                    else:
                        g.write('n%d -> n%d;\n' % (s, d))
    g.write('// solution lightpaths\n')
    # draw lightpaths
    for s in range(NET_SIZE):
        for d in range(NET_SIZE):
            if s != d:
                for w in range(WAVELENGTHS):
                    for i in range(NET_SIZE):
                        for j in range(NET_SIZE):
                            if i != j:
                                for v in range(lnks[s][d][w][i][j]):
                                    if lnks[d][s][w][j][i] > 0:
                                        g.write('n%d -> n%d [color=%d,label="n%d--n%d",dir=both];\n' % (i, j, w + 1, s, d))
                                        lnks[d][s][w][j][i] -= 1
                                    else:
                                        g.write('n%d -> n%d [color=%d,label="n%d->n%d"];\n' % (i, j, w + 1, s, d))
    g.write('}\n')
print('Run graph.dot through GraphViz for visual results.')
print('Ex: neato -Teps -o graph.eps graph.dot')