<commit_before># coding=utf-8
# cython: profile=True
# cython: linetrace=True
# cython: binding=True

"""
Wrapper around the Gurobi solver
"""

from __future__ import division, print_function

import sys

from sol.utils.const import ERR_NO_GUROBI

# Deal with different package names based on the version
try:
    from gurobi import *
except (ImportError, ModuleNotFoundError) as e:
    try:
        from gurobipy import *
    except (ImportError, ModuleNotFoundError) as e:
        print(ERR_NO_GUROBI)
        raise e

import warnings

# Ignore a warning generated by gurobi + numpy combo. gurobi checks that variable
# lists are not None and this triggers an annoying FutureWarning form numpy.
# Potentially this can be fixed by not passing it numpy arrays and instead doing
# .tolist() on each variable array (??, haven't checked that yet).
# That seems like a waste of cycles/memory, though
warnings.filterwarnings('ignore', category=FutureWarning,
                        message='comparison to `None` will result.*')

import time

import cython
from numpy import ma, zeros, arange, array, ndarray, frompyfunc, log, ones, full
from six import iterkeys, next
from six.moves import range
from cpython cimport bool
from sol.utils.exceptions import SOLException, InvalidConfigException
from sol.topology.traffic cimport TrafficClass
from sol.topology.topologynx cimport Topology
from sol.utils.logger import logger
from sol.opt.varnames cimport xp, al, be, bn, bp
from sol.utils.const import *
from sol.path.paths cimport Path, PPTC

# This is used in the node_budget function.
# Figured I'd make it private a cdef inline for speed
cdef inline _one_func(x):
    return 1

# Create a numpy ufunc for getting .x attributes of gurobi varibles
_get_x = frompyfunc(lambda x: x.x if isinstance(x, Var) else 0, 1, 1)
_is_var = frompyfunc(lambda x: isinstance(x, Var), 1, 1)

# noinspection PyClassicStyleClass
cdef class OptimizationGurobi:
    """
    Represents a SOL optimization problem.
    Uses Gurobi for building and solving the model.
    """
    def __init__(self, Topology topo, PPTC all_pptc):
        # Create a gurobi model
        self.opt = Model()
        # This will disable console output
        self.opt.params.LogToConsole = 0
        # Maximize the optimization
        self.opt.ModelSense = GRB.MAXIMIZE
        # Keep track of the topology and all traffic classes
        self.topo = topo
        self._all_pptc = all_pptc

        # Should we measure the time it takes to solve the optimization
        self._do_time = True
        self._time = 0

        # Compute the set of all resources
        res_set = set()
        for node in topo.nodes():
            rl = topo.get_resources(node)
            res_set.update(rl)
        for link in topo.links():
            rl = topo.get_resources(link)
            res_set.update(rl)
        # Ensure that number of epochs matches up across all traffic classes
        if len(set([ma.compressed(x.volFlows).size for x in all_pptc.tcs()])) != 1:
            raise ValueError(ERR_EPOCH_MISMATCH)
        # Store total number of epochs
        self.num_epochs = ma.compressed(next(all_pptc.tcs()).volFlows).size
        # Maximum number of paths for a single traffic class (to calculate the number of x_* variables)
        self._max_paths = all_pptc.max_paths(all=False)
        # The array that hold all of the x_* variables
        self._xps = zeros(shape=(all_pptc.num_tcs(), self._max_paths, self.num_epochs),
                          dtype=object)
        # Create arrays for other variables as well
        self._als = zeros((all_pptc.num_tcs(), self.num_epochs), dtype=object)
        self._bps = zeros((all_pptc.num_tcs(), self._max_paths), dtype=object)
        self._bns = zeros(topo.num_nodes(), dtype=object)
        self._bes = zeros((topo.num_nodes(), topo.num_nodes()), dtype=object)
        # Create the load dictionary -- stores linear expressions that compute load per resource
        self._load_dict = {}
        # Initialize the load dict. Maps resource->node->value
        for res in res_set:
            self._load_dict[res] = {}
            for n in topo.nodes():
                self._load_dict[res][n] = None

            for l in topo.links():
                self._load_dict[res][l] = None
        # Create all of the x_* variables
        self._add_decision_vars()
        self.allocate_flow(self._all_pptc.tcs())

        logger.debug("Initialized Gurobi wrapper")

    cdef _add_decision_vars(self):
        """
        Add desicision (x_*) variables responsible for determining
        the amount of flow on each path for each TraffiClass (in each epoch).
        """
        cdef TrafficClass tc
        cdef int epoch, pi  # path index
        # for each traffic class, path, and epoch add a variable
        # use the loops so we get different objects & names. Do not use numpy fill
        for tc in self._all_pptc.tcs():
            for pi in arange(self._all_pptc.num_paths(tc)):
                for epoch in arange(self.num_epochs):
                    self._xps[tc.ID, pi, epoch] = \
                        self.opt.addVar(lb=0, ub=1, name=xp(tc, pi, epoch))
        self.opt.update()
        logger.debug("Added desicion variables")

    cdef _add_binary_vars(self, PPTC pptc, vtypes):
        """
        Add binary types of a given type to the optimization. Supported types are variables
        per node, per link, or per path.

        :param pptc: paths per traffic class
        :param vtypes: list of variable types. See `:py:sol.utils.const.BinType`
        :return:
        """
        cdef int pi  # path index
        cdef bool mod = False  # has there been any new variables added

        # For each variable type
        for t in vtypes:
            if t == BinType.BIN_NODE:
                for n in self.topo.nodes(False):
                    if isinstance(self._bns[n], float):
                        self._bns[n] = self.opt.addVar(vtype=GRB.BINARY, name=bn(n))
                        mod = True
            elif t == BinType.BIN_EDGE:
                for u, v in self.topo.links(False):
                    if isinstance(self._bes[u, v], float):
                        self._bes[u, v] = self.opt.addVar(vtype=GRB.BINARY, name=be(u, v))
                        mod = True
            elif t == BinType.BIN_PATH:
                for tc in pptc.tcs():
                    for pi in arange(pptc.num_paths(tc)):
                        if not isinstance(self._bps[tc.ID, pi], Var):
                            self._bps[tc.ID, pi] = \
                                self.opt.addVar(vtype=GRB.BINARY, name=bp(tc, pi))
                            mod = True
            else:
                raise SOLException("Unknown binary variable type")
        # if new vars were added, update the model
        if mod:
            self.opt.update()

    # @profile
    cpdef allocate_flow(self, tcs, allocation=None):
        """
        Allocate network flow for each traffic class by allocating flow on each
        path (and summing it up for each traffic class).

        This function can be called multiple times with different traffic classes

        :param tcs: traffic classes
        :param allocation: if given, allocation for given traffic classes will
            be set to this value. Allocation must be between 0 and 1
        :raises: ValueError if the given allocation is not between 0 and 1
        """
        cdef int pi, epoch
        cdef TrafficClass tc
        logger.debug('Allocating flow')
        for tc in tcs:
            for epoch in range(self.num_epochs):
                if not isinstance(self._als[tc.ID, epoch], Var):
                    self._als[tc.ID, epoch] = v = self.opt.addVar(lb=0, ub=1, name=al(tc, epoch))
                # construct the expression: sum up all varibles per traffic class
                # self.opt.addConstr(quicksum(self._xps[tc.ID, :, epoch]) == v, name='al.{}'.format(tc.ID))
                ind = _is_var(self._xps[tc.ID, :, epoch]).astype(bool, copy=False)
                expr = LinExpr()
                expr.addTerms(ones(self._max_paths)[ind], self._xps[tc.ID, :, epoch][ind])
                self.opt.addConstr(expr, GRB.EQUAL, v)
                # If we also have an allocation value, add that constraint as well
                if allocation is not None:
                    self.opt.addConstr(v == allocation)
        # Update the model
        self.opt.update()

    cpdef route_all(self, tcs=None):
        """
        Ensure that all available traffic is routed (no drops) by forcing the allocation of flow to be 1
        for all given traffic classes.

        :param tcs: list of traffic classes
        """
        cdef TrafficClass tc
        cdef int epoch
        if tcs is None:
            tcs = self._all_pptc.tcs()
        for tc in tcs:
            for epoch in range(self.num_epochs):
                # This is sufficient as it forces both lower and upper bounds
                # on a variable to be 1
                self._als[tc.ID, epoch].lb = 1
        # Update the model
        self.opt.update()

    cpdef consume(self, tcs, unicode resource, cost_funcs, capacities, mode):
        """
        Compute the loads on a given resource by given traffic classes

        :param tcs: traffic classes that cosume the resource
        :param resource: resource to be consumed
        :param cost_funcs: an iterable containing cost functions
        :param capacities: capacities of links (or nodes)
        :param mode: either 'links' or 'nodes' -- determines what to iterate over
        """
        logger.debug(u'Consuming resource %s' % resource)
        cdef int pi  # path index
        cdef TrafficClass tc
        cdef Path path

        for tc in tcs:
            # vols = tc.volFlows.compressed()
            for pi, path in enumerate(self._all_pptc.paths(tc)):
                # print('%d %d' %(tc.ID, pi))
                if mode == NODES:  # iterate over path nodes
                    # iterate over all nodes,
                    # the cost function is responsible for determining the cost (could be 0)
                    for node in path.nodes():
                        if node in capacities and capacities[node] > 0:
                            # Make sure we don't get key error/non-existent array
                            if not isinstance(self._load_dict[resource][node], ndarray):
                                self._load_dict[resource][node] = zeros((self._all_pptc.num_tcs(),
                                                                         self._max_paths, self.num_epochs), dtype=float)
                            vals = array([func(tc, path, node) for func in cost_funcs]).max(axis=0)
                            self._load_dict[resource][node][tc.ID, pi, :] += vals / capacities[node]

                elif mode == LINKS:  # iterate over path links
                    for link in path.links():
                        if link in capacities and capacities[link] > 0:
                            if not isinstance(self._load_dict[resource][link], ndarray):
                                self._load_dict[resource][link] = zeros((self._all_pptc.num_tcs(),
                                                                         self._max_paths, self.num_epochs), dtype=float)
                            vals = array([func(tc, path, link) for func in cost_funcs]).max(axis=0)
                            self._load_dict[resource][link][tc.ID, pi, :] += vals / capacities[link]


    # @profile
    cpdef cap(self, unicode resource, caps, path_dep=False, tcs=None):
        """ Cap the usage of a given resource with a given value. 

        :param resource: the name of resource to cap
        :param tcs: traffic classes that consume the resource. If None, all traffic classes are used
        :param caps: the maximum utilization of a given resource
        """
        logger.debug("Capping resource %s" % resource)
        # Sanity check on cap value
        cdef int e, i, c
        cdef TrafficClass tc
        cdef ndarray coeffs, vars

        # if no classes are given, take all of them
        if tcs is None:
            tcs = list(self._all_pptc.tcs())
        else:
            tcs = list(tcs)
        for node_or_link in self._load_dict[resource]:
            # if no load has been computed or node/link not capped, skip
            if not isinstance(self._load_dict[resource][node_or_link], ndarray) or \
                            node_or_link not in caps:
                continue

            # For each epoch
            for e in range(self.num_epochs):
                expr = LinExpr()
                coeffs = self._load_dict[resource][node_or_link][:,:,e]
                vars = self._xps[:,:,e] if not path_dep else self._bps
                ind = coeffs.nonzero()
                # for c in range(coeffs.size):
                #     if coeffs[c] > 0:
                #         expr.add(vars[c], coeffs[c])
                expr.addTerms(coeffs[ind].reshape(-1), vars[ind].reshape(-1))
                # for tc in tcs:
                #     # extract valid x_* variables or b_* variables if resource is counted per path
                #     vars = self._xps[tc.ID, :num_paths[tc], e] if not path_dep else self._bps[tc.ID, :num_paths[tc]].reshape((1, self._max_paths, 1))
                #     # generate coefficients for variables (i.e., the load across all paths)
                #     coeffs = self._load_dict[resource][node_or_link][tc.ID, :num_paths[tc], e]
                #     ind = coeffs > 0  # get non-zero coefficients, and only add those variables
                #     expr.addTerms(coeffs[ind], vars[ind])
                self.opt.addConstr(expr, GRB.LESS_EQUAL, caps[node_or_link],
                                   # name='Cap.{}.{}.{}'.format(resource, str(node_or_link), e)
                                   )
        self.opt.update()

    cdef _req_all(self, req_type, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        # ensure x_* <= b_*
        self._disable_paths(pptc)
        cdef int pi  # path index
        cdef int u, v  # links
        # if no traffic classes, use all
        if traffic_classes is None:
            traffic_classes = pptc.keys()
        if req_type.lower() == BinType.BIN_NODE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    nodes = None
                    # depening on the mode, use either all nodes or middleboxes only
                    if node_mode == NodeConsumeMode.ALL:
                        nodes = list(path.nodes())
                    elif node_mode == NodeConsumeMode.MBOXES:
                        nodes = list(path.mboxes())
                    else:
                        raise TypeError(ERR_UNKNOWN_MODE % ('node consumption', node_mode))
                    self.opt.addConstr(self._bps[tc.ID, pi] <= self._bns[nodes])
        elif req_type.lower() == BinType.BIN_EDGE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    for u, v in path.links():
                        self.opt.addConstr(self._bps[tc, pi] <= self._bes[u, v])
        else:
            raise SOLException(u'Unknown type of constraint for _req_all()')
        # update the model
        self.opt.update()

    cdef _req_some(self, req_type, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        # ensure x_* <= b_*
        self._disable_paths(traffic_classes)
        cdef int pi  # path index
        cdef int u, v  # links
        if traffic_classes is None:
            traffic_classes = iterkeys(pptc.keys)
        if req_type.lower() == BinType.BIN_NODE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    expr = LinExpr()
                    nodes = None
                    # depening on the mode, use either all nodes or middleboxes only
                    if node_mode == NodeConsumeMode.ALL:
                        nodes = list(path.nodes())
                    elif node_mode == NodeConsumeMode.MBOXES:
                        nodes = list(path.mboxes())
                    else:
                        raise TypeError(ERR_UNKNOWN_MODE % (u'node consumption', node_mode))
                    self.opt.addConstr(self._bps[tc, pi] <= quicksum(self._bns[nodes]))
        elif req_type.lower() == BinType.BIN_EDGE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    expr = LinExpr()
                    for u, v in path.links():
                        expr.add(self._bes[u, v])
                    self.opt.addConstr(self._bps[tc, pi] <= expr)
        else:
            raise SOLException(u'Unknown type of constraint for _req_some()')
        # update the model
        self.opt.update()  # update the model

    cpdef req_all_nodes(self, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        """
        Require all nodes to be enabled for each path

        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        """
        return self._req_all(BinType.BIN_NODE, traffic_classes)

    cpdef req_all_links(self, traffic_classes=None):
        """
        Require all links to be enabled for each path

        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        """
        return self._req_all(BinType.BIN_EDGE, traffic_classes)

    cpdef req_some_nodes(self, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        """
        Require at least one node to be enabled for each path
        
        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        
        """
        return self._req_some(BinType.BIN_NODE, node_mode, traffic_classes)

    cpdef req_some_links(self, traffic_classes=None):
        """
        Require at least one link to be enabled for each path
        
        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        
        """
        return self._req_some(BinType.BIN_EDGE, traffic_classes)

    cdef _disable_paths(self, traffic_classes=None):
        """
        Add constraints which force paths where binpath_* variable is 0 to not carry any flow

        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        """
        self._add_binary_vars(self._all_pptc, [BinType.BIN_PATH])
        if traffic_classes is None:
            traffic_classes = pptc.tcs()
        cdef TrafficClass tc
        cdef Path path
        cdef int epoch
        for tc in traffic_classes:
            for pi, path in enumerate(self._all_pptc.tcs()):
                for epoch in range(self.num_epochs):
                    # x_* <= b_*
                    self.opt.addConstr(self._xps[tc.ID, pi, epoch] <= self._bps[tc.ID, pi])
        self.opt.update()

    cpdef enforce_single_path(self, traffic_classes):
        """
        Force all traffic to flow on a single path given for given traffic classes

        :param traffic_classes: traffic classes for which to enable this constraint
        """
        # ensure x_* <= b_*
        self._disable_paths(pptc)
        if traffic_classes is None:
            traffic_classes = pptc.tcs()
        cdef int pi  # path index
        for tc in traffic_classes:
            # sum of all binary path variables is 1. That is only a single 1 is allowed,
            # everything else is 0
            self.opt.addConstr(quicksum(self._bps[tc, :]) == 1)
        self.opt.update()

    cpdef flow_affinity(self, tc_pairs):
        """ 
        Ensure that the traffic for each pair of traffic classes in tc_pairs
        is processed at the same middleboxes, to ensure session-level view of the connection
        
        :param tc_pairs: list of traffic class pairs
        """

        for c1, c2 in tc_pairs:
            m2p1 = defaultdict([])
            m2p2 = defaultdict([])
            mboxes = set()
            for p, i in enumerate(self._all_pptc.paths(c1)):
                for m in p.mboxes():
                    m2p1[m].append(i)
            for p, i in enumerate(self._all_pptc.paths(c2)):
                for m in p.mboxes():
                    m2p2[m].append(i)
            for mbox in m2p1:
                for e in self.num_epochs:
                    self.opt.addConstr(quicksum(self._xps[c1.ID, m2p1[m], e]) == quicksum(self._xps[c2.ID, m2p2[m], e]))

    cpdef cap_num_paths(self, int max_paths, PPTC pptc=PPTC()):
        """
        Cap the total number of paths allowed to be enabled from the given set of paths per traffic class.
        If pptc is None, all traffic classes (and paths) are used.

        :param max_paths: number of paths
        :param pptc:
        :return:
        """
        if pptc.empty():
            pptc = self._all_pptc
        # ensure x_* <= b_*
        self._disable_paths(pptc)
        # reshapte all into a 1-d array and sum all binary path variables
        self.opt.addConstr(quicksum(self._bps.reshape(-1)) <= max_paths, name='path_cap')
        self.opt.update()  # model update

    cpdef min_latency(self, tcs=None, bool norm=True, cost_func=None, varname=None):
        """
        Add a minimize latency objective

        :param pptc: paths per traffic class
        :param norm: should the latency be normalized by the network diameter.
            Default is True, and it is better
            to keep it that way to ensure that the resulting value is in range [0, 1].
            An exception is if your cost_func already perfoms some type of normalization.
        :param cost_func: a cost function that returns the cost of a single path. If none, the *len* function
            will be used.
        :param varname: variable name to use. If None, one will be picked automatically
        :return: A gurobi model variable
        """
        if varname is None:
            varname = Objective.MIN_LATENCY.name
        if tcs is None:
            tcs = list(self._all_pptc.tcs())
        if cost_func is None:
            cost_func = len
        cdef int epoch, pi  # epoch, path index
        cdef double norm_factor = 1.0
        if norm:
            norm_factor = self.topo.diameter() * (self.topo.num_nodes() ** 2)
        per_epoch_obj = zeros(self.num_epochs, dtype=object)

        for epoch in range(self.num_epochs):
            per_epoch_obj[epoch] = latency = self.opt.addVar(name='{}_{}'.format(varname, epoch), lb=0, ub=1)
            latency_expr = LinExpr()
            for tc in tcs:
                for pi, path in enumerate(self._all_pptc.paths(tc)):
                    latency_expr.addTerms(cost_func(path)/norm_factor, self._xps[tc.ID, pi, epoch])
            # Since the global direction is maximize, and latency is minimize, we gotta do '1-' trick
            self.opt.addConstr(latency <= 1 - latency_expr)
        self.opt.update()
        return per_epoch_obj

    cpdef min_enabled_nodes(self, cost_func=None, varname=None):
        if varname is None:
            varname = Objective.MIN_ENABLED_NODES.name
        if cost_func is None:
            cost_func = _one_func
        cdef float total = 0
        cdef ndarray coeffs
        cdef int n = 0
        # if not norm:
        #     warnigns.warn(ERR_NO_NORM)
        #     total = 1
        # else:
        coeffs = array([cost_func(n) for n in range(self._bns.size)])
        total = sum(coeffs)
        for epoch in range(self.num_epochs):
            per_epoch_obj[epoch] = num_nodes = self.opt.addVar(name='{}_{}'.format(varname, epoch), lb=0, ub=1)
            self.opt.addConstr(num_nodes == 1-quicksum(self._bns * coeffs))
        self.opt.update()
        return per_epoch_obj

    cpdef node_budget(self, int bound, budget_func=None):
        """
        Enable at most *bound* nodes.

        :param bound: the maximum bound after budget_func has been computed
        :param budget_func: compute the cost of a single node.
            If None, 1 will be used, simply counting the nodes
        """
        if budget_func is None:
            budget_func = _one_func
        # double check we have binary variables
        self._add_binary_vars(None, [u'node'])
        # Sum up all binary node varibles
        expr = LinExpr()
        for n in self.topo.nodes(data=False):
            expr.addTerms(budget_func(n), self._bns[n])
        # and bound them
        self.opt.addConstr(expr <= bound)
        self.opt.update()

    # @profile
    cdef _min_load(self, unicode resource, tcs, varname):
        """
        Minimize load imposed by given traffic classes on a given resource

        :param resource: the name of resource in question
        :param tcs: the traffic classes to account for.
            If None, all traffic classes are used.
        :param varname: name of the varible. If none, one will be generated
        :return:
        """
        logger.debug('Minimizing load for %s' % resource)
        cdef int e, i
        cdef unicode per_epoch_name
        cdef int tcid, pathid
        cdef TrafficClass tc

        # sanity check traffic classes
        if tcs is None:
            tcs = list(self._all_pptc.tcs())
        else:
            tcs = list(tcs)
        # deal with variable names and create variables for each epoch
        if varname is None:
            varname = u'{}_{}'.format(MIN_LOAD_PREFIX, resource)
        per_epoch_objs = zeros(self.num_epochs, dtype=object)
        # this will create variables for objective within each epoch
        for e in range(self.num_epochs):
            per_epoch_name = u'{}_{}'.format(varname, e)
            per_epoch_objs[e] = self.opt.addVar(lb=0, ub=1, name=per_epoch_name)
        self.opt.update()  # update model

        # cdef ndarray[uint8_t, ndim=1] tcind
        cdef ndarray tcind = array([tc.ID for tc in tcs], dtype=int)
        # Proceed to compute the load
        for e in range(self.num_epochs):
            for node_or_link in self._load_dict[resource]:
                # Skip non-existent resource-node combos
                if not isinstance(self._load_dict[resource][node_or_link], ndarray):
                    continue
                expr = LinExpr()
                coeffs = self._load_dict[resource][node_or_link][tcind, :, e]
                ind = coeffs.nonzero()
                expr.addTerms(coeffs[ind].reshape(-1), self._xps[tcind, :, e][ind].reshape(-1))
                # for tc in tcs:
                #     vars = self._xps[tc.ID, :, e]
                #     # ind = _is_var(vars).astype(bool, copy=False)
                #     # get the coefficients
                #     coeffs = self._load_dict[resource][node_or_link][tc.ID, :, e]
                #     ind = coeffs > 0
                #     expr.addTerms(coeffs[ind], vars[ind])
                expr.addTerms(1, per_epoch_objs[e])
                self.opt.addConstr(expr, GRB.LESS_EQUAL, 1,
                                   # name='ml_{}_{}'.format(str(node_or_link), e)
                                   )
        self.opt.update()
        return per_epoch_objs

    cpdef min_node_load(self, unicode resource, tcs=None, varname=None):
        """
        Minimize node load for a particular resource

        :param resource: name of the resource
        :param tcs: the traffic classes that influence the objective
        :param varname: the name of the objective variable
        :return: A list of gurobi variables (one per each epoch)
        """
        return self._min_load(resource, tcs, varname)

    cpdef min_link_load(self, unicode resource, tcs=None, varname=None):
        """
        Minimize link load for a particular resource

        :param resource: name of the resource
        :param tcs: the traffic classes that influence the objective
        :param varname: now to name the objective variable. If None, a default will be provided
        :return: A list of gurobi variables (one per each epoch)
        """
        return self._min_load(resource, tcs, varname)

    cpdef max_flow(self, tcs=None, varname=None):
        """
        Maximize total network flow.

        :param tcs: traffic classes for which the flow allocation should be maximized
        :param varname: now to name the objective variable. If None, a default will be provided
        :return: A list of gurobi variables corresponding to each
        """
        # Deal with None values
        varname = Objective.MAX_FLOW.name if varname is None else varname
        if not tcs:
            tcs = list(self._all_pptc.tcs())
        # generate a variable for each epoch
        per_epoch_objs = zeros(self.num_epochs, dtype=object)
        for e in range(self.num_epochs):
            per_epoch_objs[e] = obje = self.opt.addVar(lb=0, ub=1, name='{}_{}'.format(varname, e))
            # simply sum all allocations, normalize by the number of traffic classes, and that
            # is our objective
            expr = LinExpr(full(len(tcs), 1.0/len(tcs), dtype=float), self._als[[tc.ID for tc in tcs], e])
            self.opt.addConstr(expr, GRB.GREATER_EQUAL, obje)
            # self.opt.addConstr(obje == quicksum([self._als[tc.ID, e] for tc in tcs]) / len(tcs))
        self.opt.update()
        return per_epoch_objs

    cdef _compose_obj_one_epoch(self, int epoch, ndarray obj, fairness_mode, weight_arr):
        """
        Compose multiple objectives into a single objective **for a single epoch**
        :param epoch: the epoch we are in
        :param obj: the array containing all objectives (as Gurobi vars)
        :param fairness_mode: how to compose the objectives (e.g., weighted, max-min, propfair, etc.)
        :param weight_arr: the weights array. Only applicable if the fairness_mode is "weighted"
        :return: a single gurobi variable containing the composed objectives
        """
        cdef int i
        cdef float s
        epoch_obj = self.opt.addVar(name='{}_{}'.format(THE_OBJECTIVE, epoch))
        if fairness_mode == Fairness.WEIGHTED:
            self.opt.addConstr(epoch_obj == LinExpr(weight_arr, obj))
        elif fairness_mode == Fairness.MAXMIN:
            for o in obj:
                self.opt.addConstr(epoch_obj <= o)
        elif fairness_mode == Fairness.PROPFAIR:
            # picewise appromation of proportional fairness.
            # x defines x-values at which the function should be evaluated
            x = [EPSILON, 0.01, 0.02, 0.03, 0.05, 0.08, 0.12, 0.18, 0.28, 0.43, 0.66, 1.]
            # y is the log value of x array
            y = log(x)
            # compute slopes between points
            slope = [(log(j) - log(i)) / (j - i) for i, j in zip(x, x[1:])]
            # and the y-intercepts so we can construct a line
            yint = [y[i] - slope[i] * x[i] for i in range(len(x) - 1)]
            # apply the log approximation to each objective variable individually:
            o_approx_all = []
            for o in obj:
                # Set the lower bound of o to not be 0. Because log(0) makes no sense
                # o.lb = EPSILON
                o_approx = self.opt.addVar(name='{}_{}'.format(o.VarName, 'approx'))
                # For all pieces of the linear apporoximation
                for i in range(len(x)-1):
                    s = (y[i+1] - y[i])/(x[i+1]-x[i])
                    self.opt.addConstr(o_approx <= y[i]+ o*s - s*x[i])
                # for a, b in zip(yint, slope):
                #     self.opt.addConstr(o_approx <= a + b * o)
                o_approx_all.append(o_approx)
            self.opt.addConstr(epoch_obj <= quicksum(o_approx_all))
        else:
            raise InvalidConfigException(ERR_UNKNOWN_MODE % ('fairness', fairness_mode))
        self.opt.update()
        return epoch_obj

    cpdef compose_objectives(self, ndarray obj_arr, epoch_mode, fairness_mode, weight_arr):
        """
        Compose multiple objectives, across different epochs, into a unified objective function.

        :param obj_arr: a 2-d array of gurobi variables. First dimention should be different objectives (i.e., from
            different applications). Second dimension should be objectives for a single app across epochs.
        :param weight_arr: An array containing weights of each objective function (length of this array should match
            the 1st dimension of *obj_arr*
        :param mode: the composition mode of objective functions across epochs
            (:py:class:sol.utils.const.EpochComposition). Supported values are AVG and WORST.

            AVG refers to summing up and averaging objective values across epochs, then having weights
            assigned to them.

            WORST refers to maximizing the combination of objective values for the worst epoch.
        :param weight_arr: the weights array. Only applicable if the fairness_mode is "weighted"
        :return: A gurobi variable refering to the overall objective
        """
        # TODO: Allow weight array to vary across epochs
        # TODO: stack arrays if we only have one dimension
        # if weight_arr.ndim == 1:

        the_obj = self.opt.addVar(obj=1.0, name=THE_OBJECTIVE)
        if epoch_mode == EpochComposition.AVG:
            expr = LinExpr()
            for e in range(self.num_epochs):
                # Use addTerms, it's faster
                expr.addTerms(1, self._compose_obj_one_epoch(e, obj_arr[:, e], fairness_mode, weight_arr))
            self.opt.addConstr(the_obj <= expr / self.num_epochs)
        elif epoch_mode == EpochComposition.WORST:
            for e in range(self.num_epochs):
                expr = LinExpr()
                expr.addTerms(1, self._compose_obj_one_epoch(e, obj_arr[:, e], fairness_mode, weight_arr))
                self.opt.addConstr(the_obj <= expr)
        else:
            raise ValueError(ERR_UNKNOWN_MODE % (u'epoch composition', epoch_mode))
        self.opt.update()
        return the_obj

    cpdef relax_to_lp(self):
        """
        Change integer (or binary) variables to continuous variables.

        .. warning::

            Solution is no longer guaranteed to be close to optimal (or even to make sense).
            This assumes you know what you are doing and are intentionally attempting to
            implement randomized rounding or other similar techniques.
        """
        for v in self.opt.getVars():
            if v.vType == GRB.BINARY:
                self.intvars.add(v)
                v.vType = GRB.CONTINUOUS
        self.opt.update()

    cpdef set_time_limit(self, long time):
        """
        Limit how long Gurobi looks for the solution.
        :param time: time, in milliseconds
        :return:
        """
        self.opt.params.TimeLimit = time
        self.opt.update()

    cpdef solve(self):
        """
        Solve the optimization

        ..note::
            Can be a time consuming operation
        """
        logger.debug("Running Gurobi solver")
        start = time.time()
        self.opt.optimize()
        if self._do_time:
            self._time = time.time() - start

    cpdef write(self, fname):
        """
        Writes the LP/ILP formulation to disk.
        ".lp" suffix is appended automatically

        :param fname: filename of the lp file
        """
        self.opt.write("{}.lp".format(fname))

    cpdef write_solution(self, fname):
        """
        Write the solution to disk

        :param fname: filename of the solution file.
            ".sol" suffix is appended automatically
        """
        self.opt.write("{}.sol".format(fname))

    cpdef get_gurobi_model(self):
        """
        Returns the underlying Gurobi model

        .. warning::
            Changes to the underlying model can lead to incorrect results.
            Use with caution. Advanced users only.
        """
        return self.opt

    cpdef save_hints(self, fname):
        """
        Writes the ILP hints to disk. These can be loaded later to significantly speed up solving
        of the optimization
        ".hnt" suffix is appended automatically

        :param fname: filename of the hints file
        """
        self.opt.write("{}.hnt".format(fname))

    cpdef load_hints(self, fname):
        """
        Loads previuosly saved ILP hints from disk. This can help speed up solving of the
        optimization

        :param fname: filename of the hints file
        """
        self.opt.read(self, fname)

    # cpdef get_vars(self):
    #     """
    #     Returns all Gurobi variables this optimization contains
    #     :return:
    #     """
    #     return [v.VarName for v in self.opt.getVars()]

    cpdef get_var_values(self):
        """
        Returns the mapping of variable names to values assigned by optimization
        """
        return {var.VarName: var.x for var in self.opt.getVars()}

    cpdef get_solved_objective(self, app=None):
        """
        :param app: If None, global objective is returned, otherwise the objective value for the given application
            is returned.
        :return: The objective value after the optimization is solved
        :rtype: float
        """
        if app is None:
            return self.opt.getObjective().getValue()
        else:
            return [self.opt.getVarByName('{}_{}'.format(app.name, e)).x for e in range(self.num_epochs)]

    cpdef is_solved(self):
        """
        Check if the optimization is solved

        :return: True if an optimal solution has been found
        """
        return self.opt.Status == GRB.OPTIMAL

    cpdef get_xps(self):
        """
        Return the 3-dimentional array of all decision variables.
        Dimensions are:

        1. Traffic classes (by ID)
        2. Paths
        3. Epochs


        .. warning::
            This reflects internal state of the optimization, use for read-only purposes

        :return:  numpy array
        """
        return self._xps

    cpdef get_load_dict(self):
        """
        Return the computed load for each resource and node/link (worst across all epochs)
        :return: dictionary mapping resources to nodes/link to load values
        """
        d = {}
        for resource in self._load_dict:
            d[resource] = {}
            for node_or_link in self._load_dict[resource]:
                if self._load_dict[resource][node_or_link] is None:
                    continue
                vals = _get_x(self._xps)
                loads = self._load_dict[resource][node_or_link] * vals
                # print (loads)
                d[resource][node_or_link] = loads.sum(axis=(1, 2)).max()
        return d

    cpdef get_enabled_nodes(self):
        """
        Return the list of enabled nodes, as determined by the optimization
        :return:
        """
        result = []
        for n in self.topo.nodes():
            v = self.opt.getVarByName(bn(n))
            if v is not None and v.x == 1:
                result.append(n)
        return result

    cpdef get_enabled_links(self):
        """
        Return the list of enabled links, as determined by the optimization
        :return:
        """
        result = []
        for u, v in self.topo.links():
            var = self.opt.getVarByName(be(u, v))
            if var is not None and var.x == 1:
                result.append(n)
        return result

    cpdef get_paths(self, int epoch=0):
        """
        :param epoch: Return the paths for this epoch. By default the paths for the first epoch
            are returned (that is *epoch=0*)
        :return: paths per traffic class with set flow fractions
        :rtype: :py:class:`sol.PPTC`
        """
        cdef TrafficClass tc
        cdef int pi
        cdef Path p
        cdef PPTC c = self._all_pptc.copy(deep=True)
        for tc in c.tcs():
            for pi, p in enumerate(c.paths(tc)):
                p.set_flow_fraction(self._xps[tc.ID, pi, epoch].x)
        return c


    # TODO: bring back mindiff functionality
    # TODO: add MIP hints/starts

    cpdef get_chosen_paths(self, relaxed=False):
        """
        Return the paths that were deemed "enabled" by the optimization
        :rtype: :py:class:`sol.PPTC`
        """
        # Go through the variables and pick our paths.
        cdef TrafficClass tc
        cdef int pid
        for tc in self._all_pptc.tcs():
            if relaxed:
                # Get the indices where x_p variable is 0 for all epochs, i.e., the path is unused across all epochs
                ind = [all(_get_x(x) == 0) for x in self._xps[tc.ID, :, :]]
            else:
                # Get the indices where b_p variable is 0, i.e., the path is unused.
                ind = [x.x == 0 for x in self._bps[tc.ID, :]]
            # logger.debug('-------- tc id: %d' % tc.ID)
            # logger.debug(self._all_pptc.all_paths(tc))
            # logger.debug(ind)
            # Mask the unused paths in the PPTC
            self._all_pptc.mask(tc, ind)
            # logger.debug(self._all_pptc.paths(tc))
        # print(self._xps.shape[1])
        return self._all_pptc

    cpdef fix_paths(self, PPTC pptc, fix_zero_paths=False):
        """
        Fix flow allocation of for given paths to a precise value.

        :param pptc: path per traffic class, with flow fractions set
        :param fix_zero_paths: only fix the allocation of paths where the fraction is >= 0
            (this is the default, as this function is usually used to fix a solution from a previous epoch)
        """
        cdef int pi
        for tc in pptc:
            for pi, p in enumerate(pptc[tc]):
                for e in range(self.num_epochs):
                    # Only fix non-zero paths
                    if p.flow_fraction() > 0 or not fix_zero_paths:
                        self.opt.addConstr(self._xps[tc.ID, pi, e] == \
                                           p.flow_fraction())
        self.opt.update()

    cpdef double get_time(self):
        """
        :return: The time it took to solve the optimization
        """
        return self._time

    def add_single_objective(self, name, *args, **kwargs):
        """
        Add an objective to the optimization
        :param name: name of the objective
        :param args: arguments to be passed to the objective generation function
        :param kwargs: keyword arguments to be passed to objective generation function
        :return:
        """
        epoch_objs = None
        if name == Objective.MIN_LINK_LOAD:
            epoch_objs = self.min_link_load(*args, **kwargs)
        elif name == Objective.MIN_NODE_LOAD:
            epoch_objs = self.min_node_load(*args, **kwargs)
        elif name == Objective.MIN_LATENCY:
            epoch_objs = self.min_latency(*args, **kwargs)
        elif name == Objective.MAX_FLOW:
            epoch_objs = self.max_flow(*args, **kwargs)
        elif name == Objective.MIN_ENABLED_NODES:
            epoch_objs = self.min_enabled_nodes(*args, **kwargs)
        else:
            raise InvalidConfigException("Unknown objective %s" % name)
        return epoch_objs

    def add_named_constraints(self, app):
        """
        Add supported constraints from the application to the optimization
        :param app: the application
        """
        for c in app.constraints:
            args, kwargs = c[1], c[2]
            # if c[0] == Constraint.ALLOCATE_FLOW:
            #     self.allocate_flow(*args, **kwargs)
            if c[0] == Constraint.ROUTE_ALL:
                self.route_all(*args, **kwargs)
            elif c[0] == Constraint.CAP_LINKS or cl == Constraint.CAP_NODES:
                self.cap(*args, **kwargs)
            elif c[0] == Constraint.REQ_ALL_LINKS:
                self.req_all_links(*args, **kwargs)
            elif c[0] == Constraint.REQ_ALL_NODES:
                self.req_all_nodes(*args, **kwargs)
            elif c[0] == Constraint.REQ_SOME_LINKS:
                self.req_some_links(*args, **kwargs)
            elif c[0] == Constraint.REQ_SOME_NODES:
                self.req_some_nodes(*args, **kwargs)
            elif c[0] == Constraint.FIX_PATHS:
                self.fix_paths(*args, **kwargs)
            elif c[0] == Constraint.NODE_BUDGET:
                self.node_budget(*args, **kwargs)
            else:
                raise InvalidConfigException("Unsupported constraint type %s" % c)
<commit_msg>Sped up allocate_flow() a bit by avoiding instance checks and array re-inits<commit_after># coding=utf-8
# cython: profile=True
# cython: linetrace=True
# cython: binding=True

"""
Wrapper around the Gurobi solver
"""

from __future__ import division, print_function

import sys

from sol.utils.const import ERR_NO_GUROBI

# Deal with different package names based on the version
try:
    from gurobi import *
except (ImportError, ModuleNotFoundError) as e:
    try:
        from gurobipy import *
    except (ImportError, ModuleNotFoundError) as e:
        print(ERR_NO_GUROBI)
        raise e

import warnings

# Ignore a warning generated by gurobi + numpy combo. gurobi checks that variable
# lists are not None and this triggers an annoying FutureWarning form numpy.
# Potentially this can be fixed by not passing it numpy arrays and instead doing
# .tolist() on each variable array (??, haven't checked that yet).
# That seems like a waste of cycles/memory, though
warnings.filterwarnings('ignore', category=FutureWarning,
                        message='comparison to `None` will result.*')

import time

import cython
from numpy import ma, zeros, arange, array, ndarray, frompyfunc, log, ones, full, uint8
from six import iterkeys, next
from six.moves import range
from cpython cimport bool
from sol.utils.exceptions import SOLException, InvalidConfigException
from sol.topology.traffic cimport TrafficClass
from sol.topology.topologynx cimport Topology
from sol.utils.logger import logger
from sol.opt.varnames cimport xp, al, be, bn, bp
from sol.utils.const import *
from sol.path.paths cimport Path, PPTC

# This is used in the node_budget function.
# Figured I'd make it private a cdef inline for speed
cdef inline _one_func(x):
    return 1

# Create a numpy ufunc for getting .x attributes of gurobi varibles
_get_x = frompyfunc(lambda x: x.x if isinstance(x, Var) else 0, 1, 1)
_is_var = frompyfunc(lambda x: isinstance(x, Var), 1, 1)

# noinspection PyClassicStyleClass
cdef class OptimizationGurobi:
    """
    Represents a SOL optimization problem.
    Uses Gurobi for building and solving the model.
    """
    def __init__(self, Topology topo, PPTC all_pptc):
        # Create a gurobi model
        self.opt = Model()
        # This will disable console output
        self.opt.params.LogToConsole = 0
        # Maximize the optimization
        self.opt.ModelSense = GRB.MAXIMIZE
        # Keep track of the topology and all traffic classes
        self.topo = topo
        self._all_pptc = all_pptc

        # Should we measure the time it takes to solve the optimization
        self._do_time = True
        self._time = 0

        # Compute the set of all resources
        res_set = set()
        for node in topo.nodes():
            rl = topo.get_resources(node)
            res_set.update(rl)
        for link in topo.links():
            rl = topo.get_resources(link)
            res_set.update(rl)
        # Ensure that number of epochs matches up across all traffic classes
        if len(set([ma.compressed(x.volFlows).size for x in all_pptc.tcs()])) != 1:
            raise ValueError(ERR_EPOCH_MISMATCH)
        # Store total number of epochs
        self.num_epochs = ma.compressed(next(all_pptc.tcs()).volFlows).size
        # Maximum number of paths for a single traffic class (to calculate the number of x_* variables)
        self._max_paths = all_pptc.max_paths(all=False)
        # The array that hold all of the x_* variables
        self._xps = zeros(shape=(all_pptc.num_tcs(), self._max_paths, self.num_epochs),
                          dtype=object)
        # Create arrays for other variables as well
        self._als = zeros((all_pptc.num_tcs(), self.num_epochs), dtype=object)
        self._bps = zeros((all_pptc.num_tcs(), self._max_paths), dtype=object)
        self._bns = zeros(topo.num_nodes(), dtype=object)
        self._bes = zeros((topo.num_nodes(), topo.num_nodes()), dtype=object)
        # Create the load dictionary -- stores linear expressions that compute load per resource
        self._load_dict = {}
        # Initialize the load dict. Maps resource->node->value
        for res in res_set:
            self._load_dict[res] = {}
            for n in topo.nodes():
                self._load_dict[res][n] = None

            for l in topo.links():
                self._load_dict[res][l] = None
        # Create all of the x_* variables
        self._add_decision_vars()
        self.allocate_flow(self._all_pptc.tcs())

        logger.debug("Initialized Gurobi wrapper")

    cdef _add_decision_vars(self):
        """
        Add desicision (x_*) variables responsible for determining
        the amount of flow on each path for each TraffiClass (in each epoch).
        """
        cdef TrafficClass tc
        cdef int epoch, pi  # path index
        # for each traffic class, path, and epoch add a variable
        # use the loops so we get different objects & names. Do not use numpy fill
        for tc in self._all_pptc.tcs():
            for pi in arange(self._all_pptc.num_paths(tc)):
                for epoch in arange(self.num_epochs):
                    self._xps[tc.ID, pi, epoch] = \
                        self.opt.addVar(lb=0, ub=1, name=xp(tc, pi, epoch))
        self.opt.update()
        logger.debug("Added desicion variables")

    cdef _add_binary_vars(self, PPTC pptc, vtypes):
        """
        Add binary types of a given type to the optimization. Supported types are variables
        per node, per link, or per path.

        :param pptc: paths per traffic class
        :param vtypes: list of variable types. See `:py:sol.utils.const.BinType`
        :return:
        """
        cdef int pi  # path index
        cdef bool mod = False  # has there been any new variables added

        # For each variable type
        for t in vtypes:
            if t == BinType.BIN_NODE:
                for n in self.topo.nodes(False):
                    if isinstance(self._bns[n], float):
                        self._bns[n] = self.opt.addVar(vtype=GRB.BINARY, name=bn(n))
                        mod = True
            elif t == BinType.BIN_EDGE:
                for u, v in self.topo.links(False):
                    if isinstance(self._bes[u, v], float):
                        self._bes[u, v] = self.opt.addVar(vtype=GRB.BINARY, name=be(u, v))
                        mod = True
            elif t == BinType.BIN_PATH:
                for tc in pptc.tcs():
                    for pi in arange(pptc.num_paths(tc)):
                        if not isinstance(self._bps[tc.ID, pi], Var):
                            self._bps[tc.ID, pi] = \
                                self.opt.addVar(vtype=GRB.BINARY, name=bp(tc, pi))
                            mod = True
            else:
                raise SOLException("Unknown binary variable type")
        # if new vars were added, update the model
        if mod:
            self.opt.update()

    # @profile
    cpdef allocate_flow(self, tcs, allocation=None):
        """
        Allocate network flow for each traffic class by allocating flow on each
        path (and summing it up for each traffic class).

        This function can be called multiple times with different traffic classes

        :param tcs: traffic classes
        :param allocation: if given, allocation for given traffic classes will
            be set to this value. Allocation must be between 0 and 1
        :raises: ValueError if the given allocation is not between 0 and 1
        """
        cdef int pi, epoch, np
        cdef TrafficClass tc
        cdef ndarray on
        logger.debug('Allocating flow')
        for tc in tcs:
            np = self._all_pptc.num_paths(tc)
            on = ones(np, dtype=uint8)
            for epoch in range(self.num_epochs):
                if not isinstance(self._als[tc.ID, epoch], Var):
                    self._als[tc.ID, epoch] = v = self.opt.addVar(lb=0, ub=1, name=al(tc, epoch))
                # construct the expression: sum up all varibles per traffic class
                expr = LinExpr(on, self._xps[tc.ID, :np, epoch])
                self.opt.addConstr(expr, GRB.EQUAL, v)
                # If we also have an allocation value, add that constraint as well
                if allocation is not None:
                    self.opt.addConstr(v == allocation)
        # Update the model
        self.opt.update()

    cpdef route_all(self, tcs=None):
        """
        Ensure that all available traffic is routed (no drops) by forcing the allocation of flow to be 1
        for all given traffic classes.

        :param tcs: list of traffic classes
        """
        cdef TrafficClass tc
        cdef int epoch
        if tcs is None:
            tcs = self._all_pptc.tcs()
        for tc in tcs:
            for epoch in range(self.num_epochs):
                # This is sufficient as it forces both lower and upper bounds
                # on a variable to be 1
                self._als[tc.ID, epoch].lb = 1
        # Update the model
        self.opt.update()

    cpdef consume(self, tcs, unicode resource, cost_funcs, capacities, mode):
        """
        Compute the loads on a given resource by given traffic classes

        :param tcs: traffic classes that cosume the resource
        :param resource: resource to be consumed
        :param cost_funcs: an iterable containing cost functions
        :param capacities: capacities of links (or nodes)
        :param mode: either 'links' or 'nodes' -- determines what to iterate over
        """
        logger.debug(u'Consuming resource %s' % resource)
        cdef int pi  # path index
        cdef TrafficClass tc
        cdef Path path

        for tc in tcs:
            # vols = tc.volFlows.compressed()
            for pi, path in enumerate(self._all_pptc.paths(tc)):
                # print('%d %d' %(tc.ID, pi))
                if mode == NODES:  # iterate over path nodes
                    # iterate over all nodes,
                    # the cost function is responsible for determining the cost (could be 0)
                    for node in path.nodes():
                        if node in capacities and capacities[node] > 0:
                            # Make sure we don't get key error/non-existent array
                            if not isinstance(self._load_dict[resource][node], ndarray):
                                self._load_dict[resource][node] = zeros((self._all_pptc.num_tcs(),
                                                                         self._max_paths, self.num_epochs), dtype=float)
                            vals = array([func(tc, path, node) for func in cost_funcs]).max(axis=0)
                            self._load_dict[resource][node][tc.ID, pi, :] += vals / capacities[node]

                elif mode == LINKS:  # iterate over path links
                    for link in path.links():
                        if link in capacities and capacities[link] > 0:
                            if not isinstance(self._load_dict[resource][link], ndarray):
                                self._load_dict[resource][link] = zeros((self._all_pptc.num_tcs(),
                                                                         self._max_paths, self.num_epochs), dtype=float)
                            vals = array([func(tc, path, link) for func in cost_funcs]).max(axis=0)
                            self._load_dict[resource][link][tc.ID, pi, :] += vals / capacities[link]


    # @profile
    cpdef cap(self, unicode resource, caps, path_dep=False, tcs=None):
        """ Cap the usage of a given resource with a given value. 

        :param resource: the name of resource to cap
        :param tcs: traffic classes that consume the resource. If None, all traffic classes are used
        :param caps: the maximum utilization of a given resource
        """
        logger.debug("Capping resource %s" % resource)
        # Sanity check on cap value
        cdef int e, i, c
        cdef TrafficClass tc
        cdef ndarray coeffs, vars

        # if no classes are given, take all of them
        if tcs is None:
            tcs = list(self._all_pptc.tcs())
        else:
            tcs = list(tcs)
        for node_or_link in self._load_dict[resource]:
            # if no load has been computed or node/link not capped, skip
            if not isinstance(self._load_dict[resource][node_or_link], ndarray) or \
                            node_or_link not in caps:
                continue

            # For each epoch
            for e in range(self.num_epochs):
                expr = LinExpr()
                coeffs = self._load_dict[resource][node_or_link][:,:,e]
                vars = self._xps[:,:,e] if not path_dep else self._bps
                ind = coeffs.nonzero()
                # for c in range(coeffs.size):
                #     if coeffs[c] > 0:
                #         expr.add(vars[c], coeffs[c])
                expr.addTerms(coeffs[ind].reshape(-1), vars[ind].reshape(-1))
                # for tc in tcs:
                #     # extract valid x_* variables or b_* variables if resource is counted per path
                #     vars = self._xps[tc.ID, :num_paths[tc], e] if not path_dep else self._bps[tc.ID, :num_paths[tc]].reshape((1, self._max_paths, 1))
                #     # generate coefficients for variables (i.e., the load across all paths)
                #     coeffs = self._load_dict[resource][node_or_link][tc.ID, :num_paths[tc], e]
                #     ind = coeffs > 0  # get non-zero coefficients, and only add those variables
                #     expr.addTerms(coeffs[ind], vars[ind])
                self.opt.addConstr(expr, GRB.LESS_EQUAL, caps[node_or_link],
                                   # name='Cap.{}.{}.{}'.format(resource, str(node_or_link), e)
                                   )
        self.opt.update()

    cdef _req_all(self, req_type, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        # ensure x_* <= b_*
        self._disable_paths(pptc)
        cdef int pi  # path index
        cdef int u, v  # links
        # if no traffic classes, use all
        if traffic_classes is None:
            traffic_classes = pptc.keys()
        if req_type.lower() == BinType.BIN_NODE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    nodes = None
                    # depening on the mode, use either all nodes or middleboxes only
                    if node_mode == NodeConsumeMode.ALL:
                        nodes = list(path.nodes())
                    elif node_mode == NodeConsumeMode.MBOXES:
                        nodes = list(path.mboxes())
                    else:
                        raise TypeError(ERR_UNKNOWN_MODE % ('node consumption', node_mode))
                    self.opt.addConstr(self._bps[tc.ID, pi] <= self._bns[nodes])
        elif req_type.lower() == BinType.BIN_EDGE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    for u, v in path.links():
                        self.opt.addConstr(self._bps[tc, pi] <= self._bes[u, v])
        else:
            raise SOLException(u'Unknown type of constraint for _req_all()')
        # update the model
        self.opt.update()

    cdef _req_some(self, req_type, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        # ensure x_* <= b_*
        self._disable_paths(traffic_classes)
        cdef int pi  # path index
        cdef int u, v  # links
        if traffic_classes is None:
            traffic_classes = iterkeys(pptc.keys)
        if req_type.lower() == BinType.BIN_NODE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    expr = LinExpr()
                    nodes = None
                    # depening on the mode, use either all nodes or middleboxes only
                    if node_mode == NodeConsumeMode.ALL:
                        nodes = list(path.nodes())
                    elif node_mode == NodeConsumeMode.MBOXES:
                        nodes = list(path.mboxes())
                    else:
                        raise TypeError(ERR_UNKNOWN_MODE % (u'node consumption', node_mode))
                    self.opt.addConstr(self._bps[tc, pi] <= quicksum(self._bns[nodes]))
        elif req_type.lower() == BinType.BIN_EDGE:
            for tc in traffic_classes:
                for pi, path in enumerate(pptc[tc]):
                    expr = LinExpr()
                    for u, v in path.links():
                        expr.add(self._bes[u, v])
                    self.opt.addConstr(self._bps[tc, pi] <= expr)
        else:
            raise SOLException(u'Unknown type of constraint for _req_some()')
        # update the model
        self.opt.update()  # update the model

    cpdef req_all_nodes(self, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        """
        Require all nodes to be enabled for each path

        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        """
        return self._req_all(BinType.BIN_NODE, traffic_classes)

    cpdef req_all_links(self, traffic_classes=None):
        """
        Require all links to be enabled for each path

        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        """
        return self._req_all(BinType.BIN_EDGE, traffic_classes)

    cpdef req_some_nodes(self, traffic_classes=None, node_mode=NodeConsumeMode.ALL):
        """
        Require at least one node to be enabled for each path
        
        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        
        """
        return self._req_some(BinType.BIN_NODE, node_mode, traffic_classes)

    cpdef req_some_links(self, traffic_classes=None):
        """
        Require at least one link to be enabled for each path
        
        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        :param node_mode: whether to count all nodes as enabled (:py:attr:NodeConsumeMode.ALL) or only middleboxes
            (:py:attr:`~sol.const.NodeConsumeMode.MBOXES`)
        
        """
        return self._req_some(BinType.BIN_EDGE, traffic_classes)

    cdef _disable_paths(self, traffic_classes=None):
        """
        Add constraints which force paths where binpath_* variable is 0 to not carry any flow

        :param traffic_classes: traffic classes for which to enable this constraint. If None, all are used
        """
        self._add_binary_vars(self._all_pptc, [BinType.BIN_PATH])
        if traffic_classes is None:
            traffic_classes = pptc.tcs()
        cdef TrafficClass tc
        cdef Path path
        cdef int epoch
        for tc in traffic_classes:
            for pi, path in enumerate(self._all_pptc.tcs()):
                for epoch in range(self.num_epochs):
                    # x_* <= b_*
                    self.opt.addConstr(self._xps[tc.ID, pi, epoch] <= self._bps[tc.ID, pi])
        self.opt.update()

    cpdef enforce_single_path(self, traffic_classes):
        """
        Force all traffic to flow on a single path given for given traffic classes

        :param traffic_classes: traffic classes for which to enable this constraint
        """
        # ensure x_* <= b_*
        self._disable_paths(pptc)
        if traffic_classes is None:
            traffic_classes = pptc.tcs()
        cdef int pi  # path index
        for tc in traffic_classes:
            # sum of all binary path variables is 1. That is only a single 1 is allowed,
            # everything else is 0
            self.opt.addConstr(quicksum(self._bps[tc, :]) == 1)
        self.opt.update()

    cpdef flow_affinity(self, tc_pairs):
        """ 
        Ensure that the traffic for each pair of traffic classes in tc_pairs
        is processed at the same middleboxes, to ensure session-level view of the connection
        
        :param tc_pairs: list of traffic class pairs
        """

        for c1, c2 in tc_pairs:
            m2p1 = defaultdict([])
            m2p2 = defaultdict([])
            mboxes = set()
            for p, i in enumerate(self._all_pptc.paths(c1)):
                for m in p.mboxes():
                    m2p1[m].append(i)
            for p, i in enumerate(self._all_pptc.paths(c2)):
                for m in p.mboxes():
                    m2p2[m].append(i)
            for mbox in m2p1:
                for e in self.num_epochs:
                    self.opt.addConstr(quicksum(self._xps[c1.ID, m2p1[m], e]) == quicksum(self._xps[c2.ID, m2p2[m], e]))

    cpdef cap_num_paths(self, int max_paths, PPTC pptc=PPTC()):
        """
        Cap the total number of paths allowed to be enabled from the given set of paths per traffic class.
        If pptc is None, all traffic classes (and paths) are used.

        :param max_paths: number of paths
        :param pptc:
        :return:
        """
        if pptc.empty():
            pptc = self._all_pptc
        # ensure x_* <= b_*
        self._disable_paths(pptc)
        # reshapte all into a 1-d array and sum all binary path variables
        self.opt.addConstr(quicksum(self._bps.reshape(-1)) <= max_paths, name='path_cap')
        self.opt.update()  # model update

    cpdef min_latency(self, tcs=None, bool norm=True, cost_func=None, varname=None):
        """
        Add a minimize latency objective

        :param pptc: paths per traffic class
        :param norm: should the latency be normalized by the network diameter.
            Default is True, and it is better
            to keep it that way to ensure that the resulting value is in range [0, 1].
            An exception is if your cost_func already perfoms some type of normalization.
        :param cost_func: a cost function that returns the cost of a single path. If none, the *len* function
            will be used.
        :param varname: variable name to use. If None, one will be picked automatically
        :return: A gurobi model variable
        """
        if varname is None:
            varname = Objective.MIN_LATENCY.name
        if tcs is None:
            tcs = list(self._all_pptc.tcs())
        if cost_func is None:
            cost_func = len
        cdef int epoch, pi  # epoch, path index
        cdef double norm_factor = 1.0
        if norm:
            norm_factor = self.topo.diameter() * (self.topo.num_nodes() ** 2)
        per_epoch_obj = zeros(self.num_epochs, dtype=object)

        for epoch in range(self.num_epochs):
            per_epoch_obj[epoch] = latency = self.opt.addVar(name='{}_{}'.format(varname, epoch), lb=0, ub=1)
            latency_expr = LinExpr()
            for tc in tcs:
                for pi, path in enumerate(self._all_pptc.paths(tc)):
                    latency_expr.addTerms(cost_func(path)/norm_factor, self._xps[tc.ID, pi, epoch])
            # Since the global direction is maximize, and latency is minimize, we gotta do '1-' trick
            self.opt.addConstr(latency <= 1 - latency_expr)
        self.opt.update()
        return per_epoch_obj

    cpdef min_enabled_nodes(self, cost_func=None, varname=None):
        if varname is None:
            varname = Objective.MIN_ENABLED_NODES.name
        if cost_func is None:
            cost_func = _one_func
        cdef float total = 0
        cdef ndarray coeffs
        cdef int n = 0
        # if not norm:
        #     warnigns.warn(ERR_NO_NORM)
        #     total = 1
        # else:
        coeffs = array([cost_func(n) for n in range(self._bns.size)])
        total = sum(coeffs)
        for epoch in range(self.num_epochs):
            per_epoch_obj[epoch] = num_nodes = self.opt.addVar(name='{}_{}'.format(varname, epoch), lb=0, ub=1)
            self.opt.addConstr(num_nodes == 1-quicksum(self._bns * coeffs))
        self.opt.update()
        return per_epoch_obj

    cpdef node_budget(self, int bound, budget_func=None):
        """
        Enable at most *bound* nodes.

        :param bound: the maximum bound after budget_func has been computed
        :param budget_func: compute the cost of a single node.
            If None, 1 will be used, simply counting the nodes
        """
        if budget_func is None:
            budget_func = _one_func
        # double check we have binary variables
        self._add_binary_vars(None, [u'node'])
        # Sum up all binary node varibles
        expr = LinExpr()
        for n in self.topo.nodes(data=False):
            expr.addTerms(budget_func(n), self._bns[n])
        # and bound them
        self.opt.addConstr(expr <= bound)
        self.opt.update()

    # @profile
    cdef _min_load(self, unicode resource, tcs, varname):
        """
        Minimize load imposed by given traffic classes on a given resource

        :param resource: the name of resource in question
        :param tcs: the traffic classes to account for.
            If None, all traffic classes are used.
        :param varname: name of the varible. If none, one will be generated
        :return:
        """
        logger.debug('Minimizing load for %s' % resource)
        cdef int e, i
        cdef unicode per_epoch_name
        cdef int tcid, pathid
        cdef TrafficClass tc

        # sanity check traffic classes
        if tcs is None:
            tcs = list(self._all_pptc.tcs())
        else:
            tcs = list(tcs)
        # deal with variable names and create variables for each epoch
        if varname is None:
            varname = u'{}_{}'.format(MIN_LOAD_PREFIX, resource)
        per_epoch_objs = zeros(self.num_epochs, dtype=object)
        # this will create variables for objective within each epoch
        for e in range(self.num_epochs):
            per_epoch_name = u'{}_{}'.format(varname, e)
            per_epoch_objs[e] = self.opt.addVar(lb=0, ub=1, name=per_epoch_name)
        self.opt.update()  # update model

        # cdef ndarray[uint8_t, ndim=1] tcind
        cdef ndarray tcind = array([tc.ID for tc in tcs], dtype=int)
        # Proceed to compute the load
        for e in range(self.num_epochs):
            for node_or_link in self._load_dict[resource]:
                # Skip non-existent resource-node combos
                if not isinstance(self._load_dict[resource][node_or_link], ndarray):
                    continue
                expr = LinExpr()
                coeffs = self._load_dict[resource][node_or_link][tcind, :, e]
                ind = coeffs.nonzero()
                expr.addTerms(coeffs[ind].reshape(-1), self._xps[tcind, :, e][ind].reshape(-1))
                # for tc in tcs:
                #     vars = self._xps[tc.ID, :, e]
                #     # ind = _is_var(vars).astype(bool, copy=False)
                #     # get the coefficients
                #     coeffs = self._load_dict[resource][node_or_link][tc.ID, :, e]
                #     ind = coeffs > 0
                #     expr.addTerms(coeffs[ind], vars[ind])
                expr.addTerms(1, per_epoch_objs[e])
                self.opt.addConstr(expr, GRB.LESS_EQUAL, 1,
                                   # name='ml_{}_{}'.format(str(node_or_link), e)
                                   )
        self.opt.update()
        return per_epoch_objs

    cpdef min_node_load(self, unicode resource, tcs=None, varname=None):
        """
        Minimize node load for a particular resource

        :param resource: name of the resource
        :param tcs: the traffic classes that influence the objective
        :param varname: the name of the objective variable
        :return: A list of gurobi variables (one per each epoch)
        """
        return self._min_load(resource, tcs, varname)

    cpdef min_link_load(self, unicode resource, tcs=None, varname=None):
        """
        Minimize link load for a particular resource

        :param resource: name of the resource
        :param tcs: the traffic classes that influence the objective
        :param varname: now to name the objective variable. If None, a default will be provided
        :return: A list of gurobi variables (one per each epoch)
        """
        return self._min_load(resource, tcs, varname)

    cpdef max_flow(self, tcs=None, varname=None):
        """
        Maximize total network flow.

        :param tcs: traffic classes for which the flow allocation should be maximized
        :param varname: now to name the objective variable. If None, a default will be provided
        :return: A list of gurobi variables corresponding to each
        """
        # Deal with None values
        varname = Objective.MAX_FLOW.name if varname is None else varname
        if not tcs:
            tcs = list(self._all_pptc.tcs())
        # generate a variable for each epoch
        per_epoch_objs = zeros(self.num_epochs, dtype=object)
        for e in range(self.num_epochs):
            per_epoch_objs[e] = obje = self.opt.addVar(lb=0, ub=1, name='{}_{}'.format(varname, e))
            # simply sum all allocations, normalize by the number of traffic classes, and that
            # is our objective
            expr = LinExpr(full(len(tcs), 1.0/len(tcs), dtype=float), self._als[[tc.ID for tc in tcs], e])
            self.opt.addConstr(expr, GRB.GREATER_EQUAL, obje)
            # self.opt.addConstr(obje == quicksum([self._als[tc.ID, e] for tc in tcs]) / len(tcs))
        self.opt.update()
        return per_epoch_objs

    cdef _compose_obj_one_epoch(self, int epoch, ndarray obj, fairness_mode, weight_arr):
        """
        Compose multiple objectives into a single objective **for a single epoch**
        :param epoch: the epoch we are in
        :param obj: the array containing all objectives (as Gurobi vars)
        :param fairness_mode: how to compose the objectives (e.g., weighted, max-min, propfair, etc.)
        :param weight_arr: the weights array. Only applicable if the fairness_mode is "weighted"
        :return: a single gurobi variable containing the composed objectives
        """
        cdef int i
        cdef float s
        epoch_obj = self.opt.addVar(name='{}_{}'.format(THE_OBJECTIVE, epoch))
        if fairness_mode == Fairness.WEIGHTED:
            self.opt.addConstr(epoch_obj == LinExpr(weight_arr, obj))
        elif fairness_mode == Fairness.MAXMIN:
            for o in obj:
                self.opt.addConstr(epoch_obj <= o)
        elif fairness_mode == Fairness.PROPFAIR:
            # picewise appromation of proportional fairness.
            # x defines x-values at which the function should be evaluated
            x = [EPSILON, 0.01, 0.02, 0.03, 0.05, 0.08, 0.12, 0.18, 0.28, 0.43, 0.66, 1.]
            # y is the log value of x array
            y = log(x)
            # compute slopes between points
            slope = [(log(j) - log(i)) / (j - i) for i, j in zip(x, x[1:])]
            # and the y-intercepts so we can construct a line
            yint = [y[i] - slope[i] * x[i] for i in range(len(x) - 1)]
            # apply the log approximation to each objective variable individually:
            o_approx_all = []
            for o in obj:
                # Set the lower bound of o to not be 0. Because log(0) makes no sense
                # o.lb = EPSILON
                o_approx = self.opt.addVar(name='{}_{}'.format(o.VarName, 'approx'))
                # For all pieces of the linear apporoximation
                for i in range(len(x)-1):
                    s = (y[i+1] - y[i])/(x[i+1]-x[i])
                    self.opt.addConstr(o_approx <= y[i]+ o*s - s*x[i])
                # for a, b in zip(yint, slope):
                #     self.opt.addConstr(o_approx <= a + b * o)
                o_approx_all.append(o_approx)
            self.opt.addConstr(epoch_obj <= quicksum(o_approx_all))
        else:
            raise InvalidConfigException(ERR_UNKNOWN_MODE % ('fairness', fairness_mode))
        self.opt.update()
        return epoch_obj

    cpdef compose_objectives(self, ndarray obj_arr, epoch_mode, fairness_mode, weight_arr):
        """
        Compose multiple objectives, across different epochs, into a unified objective function.

        :param obj_arr: a 2-d array of gurobi variables. First dimention should be different objectives (i.e., from
            different applications). Second dimension should be objectives for a single app across epochs.
        :param weight_arr: An array containing weights of each objective function (length of this array should match
            the 1st dimension of *obj_arr*
        :param mode: the composition mode of objective functions across epochs
            (:py:class:sol.utils.const.EpochComposition). Supported values are AVG and WORST.

            AVG refers to summing up and averaging objective values across epochs, then having weights
            assigned to them.

            WORST refers to maximizing the combination of objective values for the worst epoch.
        :param weight_arr: the weights array. Only applicable if the fairness_mode is "weighted"
        :return: A gurobi variable refering to the overall objective
        """
        # TODO: Allow weight array to vary across epochs
        # TODO: stack arrays if we only have one dimension
        # if weight_arr.ndim == 1:

        the_obj = self.opt.addVar(obj=1.0, name=THE_OBJECTIVE)
        if epoch_mode == EpochComposition.AVG:
            expr = LinExpr()
            for e in range(self.num_epochs):
                # Use addTerms, it's faster
                expr.addTerms(1, self._compose_obj_one_epoch(e, obj_arr[:, e], fairness_mode, weight_arr))
            self.opt.addConstr(the_obj <= expr / self.num_epochs)
        elif epoch_mode == EpochComposition.WORST:
            for e in range(self.num_epochs):
                expr = LinExpr()
                expr.addTerms(1, self._compose_obj_one_epoch(e, obj_arr[:, e], fairness_mode, weight_arr))
                self.opt.addConstr(the_obj <= expr)
        else:
            raise ValueError(ERR_UNKNOWN_MODE % (u'epoch composition', epoch_mode))
        self.opt.update()
        return the_obj

    cpdef relax_to_lp(self):
        """
        Change integer (or binary) variables to continuous variables.

        .. warning::

            Solution is no longer guaranteed to be close to optimal (or even to make sense).
            This assumes you know what you are doing and are intentionally attempting to
            implement randomized rounding or other similar techniques.
        """
        for v in self.opt.getVars():
            if v.vType == GRB.BINARY:
                self.intvars.add(v)
                v.vType = GRB.CONTINUOUS
        self.opt.update()

    cpdef set_time_limit(self, long time):
        """
        Limit how long Gurobi looks for the solution.
        :param time: time, in milliseconds
        :return:
        """
        self.opt.params.TimeLimit = time
        self.opt.update()

    cpdef solve(self):
        """
        Solve the optimization

        ..note::
            Can be a time consuming operation
        """
        logger.debug("Running Gurobi solver")
        start = time.time()
        self.opt.optimize()
        if self._do_time:
            self._time = time.time() - start

    cpdef write(self, fname):
        """
        Writes the LP/ILP formulation to disk.
        ".lp" suffix is appended automatically

        :param fname: filename of the lp file
        """
        self.opt.write("{}.lp".format(fname))

    cpdef write_solution(self, fname):
        """
        Write the solution to disk

        :param fname: filename of the solution file.
            ".sol" suffix is appended automatically
        """
        self.opt.write("{}.sol".format(fname))

    cpdef get_gurobi_model(self):
        """
        Returns the underlying Gurobi model

        .. warning::
            Changes to the underlying model can lead to incorrect results.
            Use with caution. Advanced users only.
        """
        return self.opt

    cpdef save_hints(self, fname):
        """
        Writes the ILP hints to disk. These can be loaded later to significantly speed up solving
        of the optimization
        ".hnt" suffix is appended automatically

        :param fname: filename of the hints file
        """
        self.opt.write("{}.hnt".format(fname))

    cpdef load_hints(self, fname):
        """
        Loads previuosly saved ILP hints from disk. This can help speed up solving of the
        optimization

        :param fname: filename of the hints file
        """
        self.opt.read(self, fname)

    # cpdef get_vars(self):
    #     """
    #     Returns all Gurobi variables this optimization contains
    #     :return:
    #     """
    #     return [v.VarName for v in self.opt.getVars()]

    cpdef get_var_values(self):
        """
        Returns the mapping of variable names to values assigned by optimization
        """
        return {var.VarName: var.x for var in self.opt.getVars()}

    cpdef get_solved_objective(self, app=None):
        """
        :param app: If None, global objective is returned, otherwise the objective value for the given application
            is returned.
        :return: The objective value after the optimization is solved
        :rtype: float
        """
        if app is None:
            return self.opt.getObjective().getValue()
        else:
            return [self.opt.getVarByName('{}_{}'.format(app.name, e)).x for e in range(self.num_epochs)]

    cpdef is_solved(self):
        """
        Check if the optimization is solved

        :return: True if an optimal solution has been found
        """
        return self.opt.Status == GRB.OPTIMAL

    cpdef get_xps(self):
        """
        Return the 3-dimentional array of all decision variables.
        Dimensions are:

        1. Traffic classes (by ID)
        2. Paths
        3. Epochs


        .. warning::
            This reflects internal state of the optimization, use for read-only purposes

        :return:  numpy array
        """
        return self._xps

    cpdef get_load_dict(self):
        """
        Return the computed load for each resource and node/link (worst across all epochs)
        :return: dictionary mapping resources to nodes/link to load values
        """
        d = {}
        for resource in self._load_dict:
            d[resource] = {}
            for node_or_link in self._load_dict[resource]:
                if self._load_dict[resource][node_or_link] is None:
                    continue
                vals = _get_x(self._xps)
                loads = self._load_dict[resource][node_or_link] * vals
                # print (loads)
                d[resource][node_or_link] = loads.sum(axis=(1, 2)).max()
        return d

    cpdef get_enabled_nodes(self):
        """
        Return the list of enabled nodes, as determined by the optimization
        :return:
        """
        result = []
        for n in self.topo.nodes():
            v = self.opt.getVarByName(bn(n))
            if v is not None and v.x == 1:
                result.append(n)
        return result

    cpdef get_enabled_links(self):
        """
        Return the list of enabled links, as determined by the optimization
        :return:
        """
        result = []
        for u, v in self.topo.links():
            var = self.opt.getVarByName(be(u, v))
            if var is not None and var.x == 1:
                result.append(n)
        return result

    cpdef get_paths(self, int epoch=0):
        """
        :param epoch: Return the paths for this epoch. By default the paths for the first epoch
            are returned (that is *epoch=0*)
        :return: paths per traffic class with set flow fractions
        :rtype: :py:class:`sol.PPTC`
        """
        cdef TrafficClass tc
        cdef int pi
        cdef Path p
        cdef PPTC c = self._all_pptc.copy(deep=True)
        for tc in c.tcs():
            for pi, p in enumerate(c.paths(tc)):
                p.set_flow_fraction(self._xps[tc.ID, pi, epoch].x)
        return c


    # TODO: bring back mindiff functionality
    # TODO: add MIP hints/starts

    cpdef get_chosen_paths(self, relaxed=False):
        """
        Return the paths that were deemed "enabled" by the optimization
        :rtype: :py:class:`sol.PPTC`
        """
        # Go through the variables and pick our paths.
        cdef TrafficClass tc
        cdef int pid
        for tc in self._all_pptc.tcs():
            if relaxed:
                # Get the indices where x_p variable is 0 for all epochs, i.e., the path is unused across all epochs
                ind = [all(_get_x(x) == 0) for x in self._xps[tc.ID, :, :]]
            else:
                # Get the indices where b_p variable is 0, i.e., the path is unused.
                ind = [x.x == 0 for x in self._bps[tc.ID, :]]
            # logger.debug('-------- tc id: %d' % tc.ID)
            # logger.debug(self._all_pptc.all_paths(tc))
            # logger.debug(ind)
            # Mask the unused paths in the PPTC
            self._all_pptc.mask(tc, ind)
            # logger.debug(self._all_pptc.paths(tc))
        # print(self._xps.shape[1])
        return self._all_pptc

    cpdef fix_paths(self, PPTC pptc, fix_zero_paths=False):
        """
        Fix flow allocation of for given paths to a precise value.

        :param pptc: path per traffic class, with flow fractions set
        :param fix_zero_paths: only fix the allocation of paths where the fraction is >= 0
            (this is the default, as this function is usually used to fix a solution from a previous epoch)
        """
        cdef int pi
        for tc in pptc:
            for pi, p in enumerate(pptc[tc]):
                for e in range(self.num_epochs):
                    # Only fix non-zero paths
                    if p.flow_fraction() > 0 or not fix_zero_paths:
                        self.opt.addConstr(self._xps[tc.ID, pi, e] == \
                                           p.flow_fraction())
        self.opt.update()

    cpdef double get_time(self):
        """
        :return: The time it took to solve the optimization
        """
        return self._time

    def add_single_objective(self, name, *args, **kwargs):
        """
        Add an objective to the optimization
        :param name: name of the objective
        :param args: arguments to be passed to the objective generation function
        :param kwargs: keyword arguments to be passed to objective generation function
        :return:
        """
        epoch_objs = None
        if name == Objective.MIN_LINK_LOAD:
            epoch_objs = self.min_link_load(*args, **kwargs)
        elif name == Objective.MIN_NODE_LOAD:
            epoch_objs = self.min_node_load(*args, **kwargs)
        elif name == Objective.MIN_LATENCY:
            epoch_objs = self.min_latency(*args, **kwargs)
        elif name == Objective.MAX_FLOW:
            epoch_objs = self.max_flow(*args, **kwargs)
        elif name == Objective.MIN_ENABLED_NODES:
            epoch_objs = self.min_enabled_nodes(*args, **kwargs)
        else:
            raise InvalidConfigException("Unknown objective %s" % name)
        return epoch_objs

    def add_named_constraints(self, app):
        """
        Add supported constraints from the application to the optimization
        :param app: the application
        """
        for c in app.constraints:
            args, kwargs = c[1], c[2]
            # if c[0] == Constraint.ALLOCATE_FLOW:
            #     self.allocate_flow(*args, **kwargs)
            if c[0] == Constraint.ROUTE_ALL:
                self.route_all(*args, **kwargs)
            elif c[0] == Constraint.CAP_LINKS or cl == Constraint.CAP_NODES:
                self.cap(*args, **kwargs)
            elif c[0] == Constraint.REQ_ALL_LINKS:
                self.req_all_links(*args, **kwargs)
            elif c[0] == Constraint.REQ_ALL_NODES:
                self.req_all_nodes(*args, **kwargs)
            elif c[0] == Constraint.REQ_SOME_LINKS:
                self.req_some_links(*args, **kwargs)
            elif c[0] == Constraint.REQ_SOME_NODES:
                self.req_some_nodes(*args, **kwargs)
            elif c[0] == Constraint.FIX_PATHS:
                self.fix_paths(*args, **kwargs)
            elif c[0] == Constraint.NODE_BUDGET:
                self.node_budget(*args, **kwargs)
            else:
                raise InvalidConfigException("Unsupported constraint type %s" % c)

{'fraction_non_alphanumeric': 0.0900470783716422, 'fraction_numerical': 0.0035890335087233453, 'mean_word_length': 1.9364733435253554, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4797771', 'n_tokens_mistral': 26369, 'n_tokens_neox': 25495, 'n_words': 8235}