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