001/*******************************************************************************
002 * Copyright (c) 2017 Pablo Pavon Marino and others.
003 * All rights reserved. This program and the accompanying materials
004 * are made available under the terms of the 2-clause BSD License 
005 * which accompanies this distribution, and is available at
006 * https://opensource.org/licenses/BSD-2-Clause
007 *
008 * Contributors:
009 *     Pablo Pavon Marino and others - initial API and implementation
010 *******************************************************************************/
011package com.net2plan.examples.ocnbook.onlineSim;
012
013
014import cern.colt.matrix.tdouble.DoubleFactory1D;
015import cern.colt.matrix.tdouble.DoubleFactory2D;
016import cern.colt.matrix.tdouble.DoubleMatrix1D;
017import cern.colt.matrix.tdouble.DoubleMatrix2D;
018import cern.jet.math.tdouble.DoublePlusMultFirst;
019import cern.jet.math.tdouble.DoublePlusMultSecond;
020import com.jom.DoubleMatrixND;
021import com.jom.OptimizationProblem;
022import com.net2plan.examples.ocnbook.offline.Offline_fa_xteFormulations;
023import com.net2plan.interfaces.networkDesign.*;
024import com.net2plan.interfaces.simulation.IEventProcessor;
025import com.net2plan.interfaces.simulation.SimEvent;
026import com.net2plan.libraries.GraphUtils.ClosedCycleRoutingException;
027import com.net2plan.utils.Constants.RoutingType;
028import com.net2plan.utils.Constants.SearchType;
029import com.net2plan.utils.*;
030
031import java.io.File;
032import java.util.*;
033
034/** 
035 * This module implements a distributed primal-decomposition-based gradient algorithm, for a coordinated adjustment of the routing in multiple domains (or cluster, or autonomous systems) in a network, so that domains do not need to exchange sensitive internal information, and minimize the average number of hops in the network.
036 *
037 * Ths event processor is adapted to permit observing the algorithm performances under user-defined conditions, 
038 * including asynchronous distributed executions, where signaling can be affected by losses and/or delays, and/or measurement errors. 
039 * The time evolution of different metrics can be stored in output files, for later processing. 
040 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec11_6_multidomainRouting_primalDecomp.m">{@code fig_sec11_6_multidomainRouting_primalDecomp.m}</a> MATLAB file used for generating the graph/s of the case study in the 
041 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm.
042 * 
043 * To simulate a network with this module, use the {@code Online_evGen_doNothing} generator.
044 * 
045 * @net2plan.keywords Multidomain network, Primal decomposition, Distributed algorithm, Flow assignment (FA), Destination-based routing, Destination-link formulation
046 * @net2plan.ocnbooksections Section 11.6
047 * @net2plan.inputParameters 
048 * @author Pablo Pavon-Marino
049 */
050public class Online_evProc_multidomainRoutingPrimalDecomp extends IEventProcessor
051{
052        private InputParameter routing_isDomainRoutingRecomputationSynchronous = new InputParameter("routing_isDomainRoutingRecomputationSynchronous", false, "true if all domains wake up synchronously to recompute its routing");
053        private InputParameter routing_averageDomainRoutingRecomputationInterval = new InputParameter("routing_averageDomainRoutingRecomputationInterval", 1.0, "Average time between two routing recomputation of a domain", 0, false, Double.MAX_VALUE, true);
054        private InputParameter routing_maxFluctuationDomainRoutingRecomputationInterval = new InputParameter("routing_maxFluctuationDomainRoutingRecomputationInterval", 0.5, "Max fluctuation in time between two routing recomputation of a domain", 0, true, Double.MAX_VALUE, true);
055        private InputParameter routing_averageMasterUpdateInterval = new InputParameter("routing_averageMasterUpdateInterval", 1.0, "Average time between two master recomputation", 0, false, Double.MAX_VALUE, true);
056        private InputParameter routing_maxFluctuationMasterUpdateInterval = new InputParameter("routing_maxFluctuationMasterUpdateInterval", 0.5, "Max fluctuation in time between two master updates", 0, true, Double.MAX_VALUE, true);
057        private InputParameter gradient_gammaStep = new InputParameter("gradient_gammaStep", 0.1, "Gamma step in the gradient algorithm", 0, false, Double.MAX_VALUE, true);
058        private InputParameter simulation_maxNumberOfMasterUpdateIntervals = new InputParameter("simulation_maxNumberOfMasterUpdateIntervals", 100.0, "Maximum number of update intervals in average per agent", 0, false, Double.MAX_VALUE, true);
059        private InputParameter simulation_randomSeed = new InputParameter("simulation_randomSeed", (long) 1, "Seed of the random number generator");
060        private InputParameter simulation_outFileNameRoot = new InputParameter("simulation_outFileNameRoot", "multidomainRoutingPrimalDecompostion", "Root of the file name to be used in the output files. If blank, no output");
061
062        private int N, E, C, Ef, D; // total number of nodes and links in the network. Ef is the number of frontier links
063        private int[] allT; // array of constants 0...N-1
064        private int[] allE; // array of constants 0...E-1
065        private double PENALIZATIONOVERSUBSCRIBED = 1E1;
066        private double PRECISIONFACTOR;
067
068        Random rng;
069
070        private Map<String, int[]> nodeGlobalIndexes_c;
071        private Map<String, int[]> frontierGlobalLinkIndexes_c;
072        private Map<String, int[]> internalAndFrontierGlobalLinkIndexes_c;
073
074        private Set<String> clusterIds;
075        private Map<Node, String> clusterIdOfNode;
076        //      private Map<Long,Integer> globalLinkIndex;
077        //      private long [] globalLinkId;
078        //      private Map<Long,Integer> globalNodeIndex;
079        //      private long [] globalNodeId;
080        private Map<String, Integer> globalClusterIndex;
081        private String[] globalClusterId;
082
083        private int[] frontierLinkGlobalIndex_ef; // one element per frontier link (in m_te and A_ce), returns its global index 
084
085        private DoubleMatrix2D global_x_te; // indexes: global node, global link 
086        private DoubleMatrix1D global_u_c; // indexes: global cluster
087        private DoubleMatrix2D global_m_te;// indexes: global node, global link (although only frontier link columns are used)
088
089        private Map<String, DoubleMatrix2D> mostUpdated_vte_fromC;
090
091        private DoubleMatrix2D global_h_nt;
092        private DoubleMatrix2D global_A_ne;
093        private DoubleMatrix1D global_u_e;
094        private DoubleMatrix2D global_h_ct;
095        private DoubleMatrix2D global_A_cef;
096
097        private Map<String, String> net2planParameters;
098
099        static final int ROUTING_MASTERUPDATE = 300;
100        static final int ROUTING_DOMAINRECOMPUTE = 301;
101
102        private TimeTrace stat_traceOf_objFun;
103        private TimeTrace stat_traceOf_x_te;
104        private TimeTrace stat_traceOf_y_e;
105
106        private NetPlan currentNetPlan;
107
108        @Override
109        public String getDescription()
110        {
111                return "This module implements a distributed primal-decomposition-based gradient algorithm, for a coordinated adjustment of the routing in multiple domains (or cluster, or autonomous systems) in a network, so that domains do not need to exchange sensitive internal information, and minimize the average number of hops in the network.";
112        }
113
114        @Override
115        public List<Triple<String, String, String>> getParameters()
116        {
117                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
118                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
119        }
120
121        @Override
122        public void initialize(NetPlan currentNetPlan, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
123        {
124                try
125                {
126                        /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
127                        InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
128
129                        if (currentNetPlan.getNumberOfLayers() != 1) throw new Net2PlanException("This algorithm works in single layer networks");
130
131                        this.rng = new Random(simulation_randomSeed.getLong());
132                        this.currentNetPlan = currentNetPlan;
133                        this.net2planParameters = net2planParameters;
134
135                        this.currentNetPlan.removeAllUnicastRoutingInformation();
136                        this.currentNetPlan.setRoutingTypeAllDemands(RoutingType.HOP_BY_HOP_ROUTING);
137
138                        this.N = this.currentNetPlan.getNumberOfNodes(); // total number of nodes in the network
139                        this.E = this.currentNetPlan.getNumberOfLinks(); // total number of nodes in the network
140                        this.D = this.currentNetPlan.getNumberOfDemands(); // total number of nodes in the network
141                        if ((E == 0) || (D == 0)) throw new Net2PlanException ("The input design should have links and demands");
142                        this.allT = new int[N];
143                        for (int n = 0; n < N; n++)
144                                allT[n] = n;
145                        this.allE = new int[E];
146                        for (int e = 0; e < E; e++)
147                                allE[e] = e;
148                        this.PRECISIONFACTOR = 1e-2;//Double.parseDouble (net2planParameters.get("precisionFactor"));
149
150                        this.clusterIds = new HashSet<String>();
151                        this.clusterIdOfNode = new HashMap<Node, String>();
152                        this.nodeGlobalIndexes_c = new HashMap<String, int[]>();
153                        this.frontierGlobalLinkIndexes_c = new HashMap<String, int[]>();
154                        this.internalAndFrontierGlobalLinkIndexes_c = new HashMap<String, int[]>();
155                        for (Node n : this.currentNetPlan.getNodes())
156                        {
157                                final int indexN = n.getIndex();
158                                final String stClusterId = n.getAttribute("clusterId");
159                                if (stClusterId == null) throw new Net2PlanException("Node " + n + " does not belong to any cluster");
160                                clusterIds.add(stClusterId);
161                                clusterIdOfNode.put(n, stClusterId);
162                                if (!nodeGlobalIndexes_c.containsKey(stClusterId))
163                                {
164                                        nodeGlobalIndexes_c.put(stClusterId, new int[] {});
165                                        frontierGlobalLinkIndexes_c.put(stClusterId, new int[] {});
166                                        internalAndFrontierGlobalLinkIndexes_c.put(stClusterId, new int[] {});
167                                }
168                                nodeGlobalIndexes_c.put(stClusterId, IntUtils.concatenate(nodeGlobalIndexes_c.get(stClusterId), new int[] { indexN }));
169                        }
170
171                        this.C = nodeGlobalIndexes_c.size();
172                        this.globalClusterIndex = new HashMap<String, Integer>();
173                        this.globalClusterId = new String[C];
174                        int counter = 0;
175                        for (String cId : clusterIds)
176                        {
177                                globalClusterId[counter] = cId;
178                                globalClusterIndex.put(cId, counter);
179                                counter++;
180                        }
181
182                        /* For each cluster, the first links in xte are the frontier links, in the same order as they appear in cluster mte */
183                        this.frontierLinkGlobalIndex_ef = new int[0];
184                        for (Link e : this.currentNetPlan.getLinks())
185                        {
186                                final int indexE = e.getIndex();
187                                final Node a_e = e.getOriginNode();
188                                final Node b_e = e.getDestinationNode();
189                                final String clusterId_ae = clusterIdOfNode.get(a_e);
190                                final String clusterId_be = clusterIdOfNode.get(b_e);
191                                if (!clusterId_ae.equals(clusterId_be))
192                                {
193                                        frontierGlobalLinkIndexes_c.put(clusterId_ae, IntUtils.concatenate(frontierGlobalLinkIndexes_c.get(clusterId_ae), new int[] { indexE }));
194                                        frontierGlobalLinkIndexes_c.put(clusterId_be, IntUtils.concatenate(frontierGlobalLinkIndexes_c.get(clusterId_be), new int[] { indexE }));
195                                        internalAndFrontierGlobalLinkIndexes_c.put(clusterId_ae, IntUtils.concatenate(internalAndFrontierGlobalLinkIndexes_c.get(clusterId_ae), new int[] { indexE }));
196                                        internalAndFrontierGlobalLinkIndexes_c.put(clusterId_be, IntUtils.concatenate(internalAndFrontierGlobalLinkIndexes_c.get(clusterId_be), new int[] { indexE }));
197                                        frontierLinkGlobalIndex_ef = IntUtils.concatenate(frontierLinkGlobalIndex_ef, new int[] { indexE });
198                                }
199                        }
200                        for (Link e : this.currentNetPlan.getLinks())
201                        {
202                                final int indexE = e.getIndex();
203                                final Node a_e = e.getOriginNode();
204                                final Node b_e = e.getDestinationNode();
205                                final String clusterId_ae = clusterIdOfNode.get(a_e);
206                                final String clusterId_be = clusterIdOfNode.get(b_e);
207                                if (clusterId_ae.equals(clusterId_be)) internalAndFrontierGlobalLinkIndexes_c.put(clusterId_ae, IntUtils.concatenate(internalAndFrontierGlobalLinkIndexes_c.get(clusterId_ae), new int[] { indexE }));
208                        }
209
210                        this.Ef = frontierLinkGlobalIndex_ef.length;
211
212                        this.global_x_te = DoubleFactory2D.sparse.make(N, E);
213                        this.global_u_c = DoubleFactory1D.dense.make(C);
214                        this.global_m_te = DoubleFactory2D.dense.make(N, E);
215                        this.mostUpdated_vte_fromC = new HashMap<String, DoubleMatrix2D>();
216                        for (String cId : clusterIds)
217                                mostUpdated_vte_fromC.put(cId, null);
218                        this.global_A_ne = this.currentNetPlan.getMatrixNodeLinkIncidence();
219                        this.global_u_e = this.currentNetPlan.getVectorLinkCapacity();
220                        this.global_h_nt = this.currentNetPlan.getMatrixNode2NodeOfferedTraffic();
221                        for (int t = 0; t < N; t++)
222                        {
223                                double h_t = 0;
224                                for (int n = 0; n < N; n++)
225                                        h_t += global_h_nt.get(n, t);
226                                global_h_nt.set(t, t, -h_t);
227                        }
228
229                        this.global_h_ct = DoubleFactory2D.dense.make(C, N);
230                        for (String clusterId : clusterIds)
231                        {
232                                final int index_c = globalClusterIndex.get(clusterId);
233                                for (Node t : this.currentNetPlan.getNodes())
234                                {
235                                        final int index_t = t.getIndex();
236                                        final int index_cluster_t = globalClusterIndex.get(clusterIdOfNode.get(t));
237                                        final boolean targetInsideTheCluster = (index_c == index_cluster_t);
238                                        if (targetInsideTheCluster)
239                                        {
240                                                /* h_ct is what other nodes outside cluster generate towards t (in negative) */
241                                                for (Node n : this.currentNetPlan.getNodes())
242                                                {
243                                                        final int index_n = n.getIndex();
244                                                        final int index_cluster_n = globalClusterIndex.get(clusterIdOfNode.get(n));
245                                                        if (index_cluster_n == index_c) continue; // only outside cluster nodes matter
246                                                        global_h_ct.set(index_c, index_t, global_h_ct.get(index_c, index_t) - global_h_nt.get(index_n, index_t));
247                                                }
248                                        } else
249                                        {
250                                                /* h_ct is what this cluster generate towards t (in negative) */
251                                                for (int index_n : nodeGlobalIndexes_c.get(clusterId))
252                                                {
253                                                        final int index_cluster_n = globalClusterIndex.get(clusterIdOfNode.get(currentNetPlan.getNode(index_n)));
254                                                        if (index_cluster_n != index_c) throw new RuntimeException("Bad");
255                                                        global_h_ct.set(index_c, index_t, global_h_ct.get(index_c, index_t) + global_h_nt.get(index_n, index_t));
256                                                }
257                                        }
258                                }
259                        }
260
261                        this.global_A_cef = DoubleFactory2D.dense.make(C, Ef);
262                        for (int index_ef = 0; index_ef < frontierLinkGlobalIndex_ef.length; index_ef++)
263                        {
264                                final int index_e = frontierLinkGlobalIndex_ef[index_ef];
265                                final Link e = currentNetPlan.getLink(index_e);
266                                final Node a_e = e.getOriginNode();
267                                final Node b_e = e.getDestinationNode();
268                                final int indexC_ae = this.globalClusterIndex.get(clusterIdOfNode.get(a_e));
269                                final int indexC_be = this.globalClusterIndex.get(clusterIdOfNode.get(b_e));
270                                if (indexC_ae != indexC_be)
271                                {
272                                        global_A_cef.set(indexC_ae, index_ef, 1.0);
273                                        global_A_cef.set(indexC_be, index_ef, -1.0);
274                                }
275                        }
276
277                        stat_traceOf_objFun = new TimeTrace();
278                        stat_traceOf_x_te = new TimeTrace();
279                        stat_traceOf_y_e = new TimeTrace();
280
281                        /* Initialize the master variables. The feasible, closest to zero */
282                        this.global_m_te = projectMaster(global_m_te, null);
283
284                        /* Create the cluster topologies */
285                        /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
286                        for (String cId : clusterIds)
287                                scheduleEvent(new SimEvent(0, SimEvent.DestinationModule.EVENT_PROCESSOR, ROUTING_DOMAINRECOMPUTE, cId));
288                        scheduleEvent(new SimEvent(0.0001, SimEvent.DestinationModule.EVENT_PROCESSOR, ROUTING_MASTERUPDATE, -1));
289
290                        checksInitialize();
291                } catch (Exception e)
292                {
293                        e.printStackTrace();
294                        throw e;
295                }
296        }
297
298        @Override
299        public void processEvent(NetPlan currentNetPlan, SimEvent event)
300        {
301                        final double time = event.getEventTime();
302                        switch (event.getEventType())
303                        {
304                        case ROUTING_DOMAINRECOMPUTE:
305                        {
306                                final String clusterId = (String) event.getEventObject();
307                                final int index_c = globalClusterIndex.get(clusterId);
308
309                                /* Takes the most updated master information it has */
310                                Triple<DoubleMatrix2D, DoubleMatrix2D, Double> p = computeSubproblem(clusterId, this.currentNetPlan);
311                                DoubleMatrix2D cluster_xte = p.getFirst();
312                                DoubleMatrix2D cluster_vte = p.getSecond();
313                                double u_c = p.getThird();
314
315                                final int[] cluster_Eif = internalAndFrontierGlobalLinkIndexes_c.get(clusterId);
316
317                                /* global x_te can be not feasible, if different domains operate asynchornously */
318                                global_x_te.viewSelection(allT, cluster_Eif).assign(cluster_xte.copy().viewSelection(allT, cluster_Eif));
319                                global_u_c.set(index_c, u_c);
320                                mostUpdated_vte_fromC.put(clusterId, cluster_vte.copy());
321
322                                /* They will be the new forwarding rules */
323                                DoubleMatrix2D new_fde = this.currentNetPlan.getMatrixDemandBasedForwardingRules();
324
325                                /* Remove all the forwarding rules of the cluster nodes */
326                                for (int index_n : nodeGlobalIndexes_c.get(clusterId))
327                                        for (Link e : currentNetPlan.getNode(index_n).getOutgoingLinks())
328                                                new_fde.viewColumn(e.getIndex()).assign(0);
329
330                                /* Create new forwarding rules according to the routing computed */
331                                for (int index_n : nodeGlobalIndexes_c.get(clusterId))
332                                        for (Node t : this.currentNetPlan.getNodes())
333                                        {
334                                                final Node n = this.currentNetPlan.getNode(index_n);
335                                                if (n == t) continue;
336                                                final int index_t = t.getIndex();
337                                                final Set<Demand> demandIdsEndingIn_t = t.getIncomingDemands();
338                                                double outTraffic_nt = 0;
339                                                for (Link e : n.getOutgoingLinks())
340                                                        outTraffic_nt += global_x_te.get(index_t, e.getIndex());
341
342                                                Link outLinkMaxFte = null;
343                                                double maxFte = -1;
344                                                double accumFteToForward = 0;
345                                                for (Link e : n.getOutgoingLinks())
346                                                {
347                                                        final double f_te = (outTraffic_nt == 0) ? 0 : global_x_te.get(index_t, e.getIndex()) / outTraffic_nt;
348                                                        if (f_te > maxFte)
349                                                        {
350                                                                maxFte = f_te;
351                                                                outLinkMaxFte = e;
352                                                        }
353                                                        if (f_te > PRECISIONFACTOR)
354                                                                for (Demand d : demandIdsEndingIn_t)
355                                                                        new_fde.set(d.getIndex(), e.getIndex(), f_te); //changingFR.put (Pair.of (d,e),f_te);
356                                                        else
357                                                                accumFteToForward += f_te;
358                                                }
359                                                for (Demand d : demandIdsEndingIn_t)
360                                                        new_fde.set(d.getIndex(), outLinkMaxFte.getIndex(), new_fde.get(d.getIndex(), outLinkMaxFte.getIndex()) + accumFteToForward);
361                                        }
362
363                                DoubleMatrix2D outFR_dn = new_fde.zMult(currentNetPlan.getMatrixNodeLinkOutgoingIncidence().viewDice(), null);
364                                for (int n = 0; n < N; n++)
365                                        for (int d = 0; d < D; d++)
366                                                if ((Math.abs(outFR_dn.get(d, n) - 1) > 1e-3) && (Math.abs(outFR_dn.get(d, n)) > 1e-3)) throw new RuntimeException("Bad. outFR_dn:  " + outFR_dn);
367
368                                /* In some rare occasions, the routing may create loops. These routings are not applied to the network */
369                                try { this.currentNetPlan.setForwardingRules(new_fde , null); } catch (ClosedCycleRoutingException e) { System.out.println ("The routing after this domain update would produce a closed loop. It is not applied.");    }
370
371                                if (time > 0)
372                                {
373                                        DoubleMatrix1D blockedTraffic = currentNetPlan.getVectorDemandBlockedTraffic();
374                                        DoubleMatrix1D linkOversubscription = currentNetPlan.getVectorLinkOversubscribedTraffic();
375                                        if (blockedTraffic.getMaxLocation()[0] > 1e-3) System.out.println ("Some traffic is lost. Sum blocked traffic of the demands: " + blockedTraffic.zSum());
376                                        if (linkOversubscription.getMaxLocation()[0] > 1e-3) System.out.println ("Some links are ovsersubscribed. Sum oversubscription traffic in the links: " + linkOversubscription.zSum());
377                                        stat_traceOf_objFun.add(time, computeObjectiveFunctionFromNetPlan(this.currentNetPlan));
378                                        stat_traceOf_x_te.add(time, this.currentNetPlan.getMatrixDestination2LinkTrafficCarried());
379                                        stat_traceOf_y_e.add(time, this.currentNetPlan.getVectorLinkCarriedTraffic());
380                                }
381
382                                if (time > this.simulation_maxNumberOfMasterUpdateIntervals.getDouble() * routing_averageMasterUpdateInterval.getDouble())
383                                {
384                                        this.endSimulation();
385                                }
386
387                                final double nextWakeUpTime = (routing_isDomainRoutingRecomputationSynchronous.getBoolean()) ? time + routing_averageDomainRoutingRecomputationInterval.getDouble() : Math.max(time, time + routing_averageDomainRoutingRecomputationInterval.getDouble() + routing_maxFluctuationDomainRoutingRecomputationInterval.getDouble() * (rng.nextDouble() - 0.5));
388                                this.scheduleEvent(new SimEvent(nextWakeUpTime, SimEvent.DestinationModule.EVENT_PROCESSOR, ROUTING_DOMAINRECOMPUTE, clusterId));
389
390                                break;
391                        }
392                        case ROUTING_MASTERUPDATE: // A node receives from an out neighbor the q_nt for any destination
393                        {
394                                DoubleMatrix2D gradient_te = DoubleFactory2D.dense.make(N, Ef, 1.0); // initially all ones
395
396                                for (String clusterId : clusterIds)
397                                {
398                                        DoubleMatrix2D v_te_fromC = mostUpdated_vte_fromC.get(clusterId);
399                                        gradient_te.assign(v_te_fromC.viewSelection(allT, frontierLinkGlobalIndex_ef), DoublePlusMultFirst.plusMult(1)); // -1 no
400                                }
401
402                                //final double gamma = gradient_gammaStep.getDouble() / (1+time*time);
403                                final double gamma = gradient_gammaStep.getDouble();
404                                //          final double gamma = gradient_gammaStep.getDouble() / (1+time);
405                                //          final double gamma = time < 20? gradient_gammaStep.getDouble() : gradient_gammaStep.getDouble() / (Math.sqrt(1+time));
406                                DoubleMatrix2D candidate_mte = global_m_te.copy();
407                                candidate_mte.viewSelection(allT, frontierLinkGlobalIndex_ef).assign(gradient_te, DoublePlusMultSecond.plusMult(-gamma));
408
409                                /* Now we have to project the candidate master variables into the set of feasible master variables */
410                                this.global_m_te = projectMaster(candidate_mte, global_m_te);
411
412                                if (time > this.simulation_maxNumberOfMasterUpdateIntervals.getDouble() * routing_averageDomainRoutingRecomputationInterval.getDouble())
413                                {
414                                        endSimulation();
415                                }
416
417                                final double nextMasterTime = Math.max(time, time + routing_averageMasterUpdateInterval.getDouble() + routing_maxFluctuationMasterUpdateInterval.getDouble() * (rng.nextDouble() - 0.5));
418                                scheduleEvent(new SimEvent(nextMasterTime, SimEvent.DestinationModule.EVENT_PROCESSOR, ROUTING_MASTERUPDATE, -1));
419                                break;
420                        }
421
422                        default:
423                                throw new RuntimeException("Unexpected received event");
424                        }
425        }
426
427        public String finish(StringBuilder st, double simTime)
428        {
429                if (simulation_outFileNameRoot.getString().equals("")) return null;
430                stat_traceOf_objFun.printToFile(new File(simulation_outFileNameRoot.getString() + "_objFunc.txt"));
431                stat_traceOf_x_te.printToFile(new File(simulation_outFileNameRoot.getString() + "_xte.txt"));
432                stat_traceOf_y_e.printToFile(new File(simulation_outFileNameRoot.getString() + "_ye.txt"));
433                final NetPlan optNetPlan = computeOptimumSolution();
434                TimeTrace.printToFile(new File(simulation_outFileNameRoot.getString() + "_jom_objFunc.txt"), new double[] { computeObjectiveFunctionFromNetPlan(optNetPlan) });
435                TimeTrace.printToFile(new File(simulation_outFileNameRoot.getString() + "_jom_xte.txt"), optNetPlan.getMatrixDestination2LinkTrafficCarried());
436                TimeTrace.printToFile(new File(simulation_outFileNameRoot.getString() + "_jom_ye.txt"), optNetPlan.getVectorLinkCarriedTraffic());
437                return null;
438        }
439
440        /* They use the most updated m_te information */
441        private Triple<DoubleMatrix2D, DoubleMatrix2D, Double> computeSubproblem(String clusterId, NetPlan np)
442        {
443                final int[] cluster_Ef = frontierGlobalLinkIndexes_c.get(clusterId);
444                final int[] cluster_N = nodeGlobalIndexes_c.get(clusterId);
445
446                DoubleMatrix2D h_nt = global_h_nt.copy().viewSelection(cluster_N, allT);
447                DoubleMatrix2D A_ne = global_A_ne.copy().viewSelection(cluster_N, allE);
448                DoubleMatrix1D u_e = global_u_e.copy();
449
450                OptimizationProblem op = new OptimizationProblem();
451
452                op.setInputParameter("h_nt", h_nt.toArray());
453                op.setInputParameter("A_ne", A_ne.toArray());
454                op.setInputParameter("u_e", u_e.toArray(), "row");
455                op.setInputParameter("m_te", global_m_te.toArray());
456                op.setInputParameter("allT", allT, "row");
457                op.setInputParameter("frontC", cluster_Ef, "row");
458                op.setInputParameter("M", PENALIZATIONOVERSUBSCRIBED);
459
460                op.addDecisionVariable("x_te", false, new int[] { N, E }, 0, Double.MAX_VALUE);
461                op.addDecisionVariable("u", false, new int[] { 1, 1 }, 0, Double.MAX_VALUE);
462
463                /* Set some input parameters */
464                op.setObjectiveFunction("minimize", "sum(x_te) + M*u");
465
466                op.addConstraint("A_ne * (x_te') == h_nt"); // the flow-conservation constraints (NxN constraints)
467                op.addConstraint("sum(x_te,1) <= u_e + u"); // the capacity constraints (E constraints)
468                op.addConstraint("x_te(allT,frontC) == m_te(allT,frontC)", "v_te"); // the capacity constraints (E constraints)
469
470                /* Call the solver to solve the problem */
471                op.solve("cplex");
472
473                /* If an optimal solution was not found, quit */
474                if (!op.solutionIsOptimal()) throw new RuntimeException("An optimal solution was not found");
475
476                /* Retrieve the optimum solutions. Convert the bps into Erlangs */
477                DoubleMatrix2D x_te = op.getPrimalSolution("x_te").view2D();
478                double u = op.getPrimalSolution("u").toValue();
479
480                /* v_te: zeros in no frontier links, the appropriate values in this cluster frontier links */
481                DoubleMatrix2D v_te = DoubleFactory2D.dense.make(N, E);
482                v_te.viewSelection(allT, cluster_Ef).assign(op.getMultipliersOfConstraint("v_te").view2D());
483
484                return Triple.of(x_te, v_te, u);
485        }
486
487        private DoubleMatrix2D projectMaster(DoubleMatrix2D can_mte, DoubleMatrix2D old_mte)
488        {
489                OptimizationProblem op = new OptimizationProblem();
490
491                op.addDecisionVariable("m_te", false, new int[] { N, Ef }, 0, Double.MAX_VALUE);
492                op.setInputParameter("h_ct", global_h_ct.toArray());
493                op.setInputParameter("A_ce", global_A_cef.toArray());
494                op.setInputParameter("can_mte", can_mte.viewSelection(allT, frontierLinkGlobalIndex_ef).toArray());
495                op.setInputParameter("u_e", global_u_e.viewSelection(frontierLinkGlobalIndex_ef).toArray(), "row");
496
497                /* Set some input parameters */
498                op.setObjectiveFunction("minimize", "sum((m_te-can_mte)^2)");
499                op.addConstraint("A_ce * (m_te') == h_ct"); // the flow-conservation constraints (NxN constraints)
500                op.addConstraint("sum(m_te,1) <= u_e"); // the capacity constraints (E constraints)
501
502                for (String cId : clusterIds)
503                {
504                        final int indexC = this.globalClusterIndex.get(cId);
505                        final int[] nIdsCluster = nodeGlobalIndexes_c.get(cId);
506                        final int[] frontierOutLinkIndexesInMte = DoubleUtils.find(global_A_cef.viewRow(indexC).toArray(), 1.0, SearchType.ALL);
507                        op.setInputParameter("destCluster", nIdsCluster, "row");
508                        op.setInputParameter("clusterOutLinks", frontierOutLinkIndexesInMte, "row");
509                        op.addConstraint("m_te(destCluster,clusterOutLinks) == 0"); // traffic to a cluster does not leave the cluster
510                }
511                if (old_mte != null) op.setInitialSolution("m_te", new DoubleMatrixND(old_mte.viewSelection(allT, frontierLinkGlobalIndex_ef).copy()));
512
513                /* Call the solver to solve the problem */
514                try
515                {
516                        op.solve("ipopt");
517                } catch (Exception e)
518                {
519                        e.printStackTrace();
520                } // sometimes IPOPT is unstable. do not stop the algorithm then
521
522                /* If an optimal solution was not found, quit */
523                if (!op.solutionIsOptimal()) return old_mte.copy(); // sometimes IPOPT is unstable. do not stop the algorithm then
524
525                /* Retrieve the optimum solutions. Convert the bps into Erlangs */
526                double[][] m_te_array = (double[][]) op.getPrimalSolution("m_te").toArray();
527
528                for (String cId : clusterIds)
529                {
530                        final int indexC = this.globalClusterIndex.get(cId);
531                        final int[] nIdsCluster = nodeGlobalIndexes_c.get(cId);
532                        final int[] frontierOutLinkIndexesInMte = DoubleUtils.find(global_A_cef.viewRow(indexC).toArray(), 1.0, SearchType.ALL);
533                        for (int t : nIdsCluster)
534                                for (int e : frontierOutLinkIndexesInMte)
535                                        if (m_te_array[t][e] > PRECISIONFACTOR)
536                                                throw new RuntimeException("Bad: m_te_array[t][e]: " + m_te_array[t][e]);
537                                        else
538                                                m_te_array[t][e] = 0;
539                }
540
541                DoubleMatrix2D m_te = DoubleFactory2D.dense.make(N, E);
542                m_te.viewSelection(allT, frontierLinkGlobalIndex_ef).assign(m_te_array);
543                return m_te;
544        }
545
546        private double computeObjectiveFunctionFromNetPlan(NetPlan np)
547        {
548                double optCost = 0;
549                double maxOverssubs = 0;
550                for (Link e : np.getLinks())
551                {
552                        final double y_e = e.getCarriedTraffic();
553                        final double u_e = e.getCapacity();
554                        maxOverssubs = Math.max(maxOverssubs, y_e - u_e);
555                        optCost += Math.min(y_e, u_e);
556                }
557                optCost += PENALIZATIONOVERSUBSCRIBED * maxOverssubs;
558                return optCost;
559        }
560
561        private NetPlan computeOptimumSolution()
562        {
563                /* Compute the optimum solution */
564                Map<String, String> m = new HashMap<String, String>();
565                m.put("optimizationTarget", "min-av-num-hops");
566                m.put("solverName", "cplex"); // use any linear solver reachable from JOM
567                m.put("solverLibraryName", "");
568                m.put("maxSolverTimeInSeconds", "-1");
569                m.put("binaryRatePerTrafficUnit_bps", "1");
570                m.put("averagePacketLengthInBytes", "1");
571                m.put("nonBifurcatedRouting", "false");
572                NetPlan npCopy = this.currentNetPlan.copy();
573                new Offline_fa_xteFormulations().executeAlgorithm(npCopy, m, this.net2planParameters);
574                return npCopy;
575        }
576
577        private void checksInitialize()
578        {
579                /* Check in h_ct, the columns sum zero */
580                for (int t = 0; t < this.currentNetPlan.getNumberOfNodes(); t++)
581                        if (Math.abs(global_h_ct.viewColumn(t).zSum()) > PRECISIONFACTOR) throw new RuntimeException("Bad");
582                /* Check in h_nt, the columns sum zero */
583                for (int t = 0; t < this.currentNetPlan.getNumberOfNodes(); t++)
584                        if (Math.abs(global_h_nt.viewColumn(t).zSum()) > PRECISIONFACTOR) throw new RuntimeException("Bad");
585        }
586
587}