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
012import java.io.File;
013import java.io.FileOutputStream;
014import java.io.PrintStream;
015import java.util.Arrays;
016import java.util.HashMap;
017import java.util.LinkedList;
018import java.util.List;
019import java.util.Map;
020import java.util.Random;
021
022import cern.colt.matrix.tdouble.DoubleFactory2D;
023import cern.colt.matrix.tdouble.DoubleMatrix1D;
024import cern.colt.matrix.tdouble.DoubleMatrix2D;
025import cern.colt.matrix.tint.IntFactory3D;
026import cern.colt.matrix.tint.IntMatrix3D;
027import cern.jet.random.tdouble.Poisson;
028
029import com.jom.OptimizationProblem;
030import com.net2plan.interfaces.networkDesign.Demand;
031import com.net2plan.interfaces.networkDesign.Link;
032import com.net2plan.interfaces.networkDesign.Net2PlanException;
033import com.net2plan.interfaces.networkDesign.NetPlan;
034import com.net2plan.interfaces.networkDesign.Node;
035import com.net2plan.interfaces.networkDesign.Route;
036import com.net2plan.interfaces.simulation.IEventProcessor;
037import com.net2plan.interfaces.simulation.SimEvent;
038import com.net2plan.libraries.GraphUtils;
039import com.net2plan.utils.Constants.RoutingType;
040import com.net2plan.utils.DoubleUtils;
041import com.net2plan.utils.InputParameter;
042import com.net2plan.utils.IntUtils;
043import com.net2plan.utils.Pair;
044import com.net2plan.utils.TimeTrace;
045import com.net2plan.utils.Triple;
046
047/** 
048 * This module implements a distributed dual-gradient based algorithm for adapting the network routing to the one which minimizes the average number of hops, that results in a purely decentralized backpressure scheme.
049 *
050 * Ths event processor is adapted to permit observing the algorithm performances under user-defined conditions, 
051 * including asynchronous distributed executions, where signaling can be affected by losses and/or delays, and/or measurement errors. 
052 * The time evolution of different metrics can be stored in output files, for later processing. 
053 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec10_3_backpressureRoutingDual.m">{@code fig_sec10_3_backpressureRoutingDual.m}</a> MATLAB file used for generating the graph/s of the case study in the 
054 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm.
055 * 
056 * To simulate a network with this module, use the {@code Online_evGen_doNothing} generator.
057 * 
058 * @net2plan.keywords Flow assignment (FA), Distributed algorithm, Backpressure routing, Dual gradient algorithm
059 * @net2plan.ocnbooksections Section 10.3
060 * @net2plan.inputParameters 
061 * @author Pablo Pavon-Marino
062 */
063public class Online_evProc_backpressureRoutingDual extends IEventProcessor 
064{
065        private static PrintStream getNulFile () { try { return new PrintStream (new FileOutputStream ("NUL") , false); } catch (Exception e) {e.printStackTrace(); throw new RuntimeException ("Not NUL file"); }   } 
066        private PrintStream log = getNulFile (); //System.err;//new PrintStream (new FileOutputStream ("NUL") , true);
067
068        private Random rng;
069        private int N,E,D;
070        
071        private static final int SIGNALING_WAKEUPTOSENDMESSAGE = 300; // A NODE SIGNAL QUEUES TO NEIGHBOR NODES
072        private static final int SIGNALING_RECEIVEDMESSAGE = 301;  // A NODE RECEIVED SIGNAL QUEUES TO NEIGHBOR NODES
073        private static final int UPDATE_INITTIMESLOTFORSCHEDULINGPACKETS = 302; // ALL NODES SCHEDULE ONE PACKET PER LINK
074        private static final int UPDATE_STATISTICTRACES = 303; // STAT UPDATE OF N2P
075
076        private InputParameter signaling_isSynchronous = new InputParameter ("signaling_isSynchronous", false , "true if all the distributed agents involved wake up synchronously to send the signaling messages");
077        private InputParameter signaling_averageInterMessageTime = new InputParameter ("signaling_averageInterMessageTime", 1.0 , "Average time between two signaling messages sent by an agent" , 0 , false , Double.MAX_VALUE , true);
078        private InputParameter signaling_maxFluctuationInterMessageTime = new InputParameter ("signaling_maxFluctuationInterMessageTime", 0.5 , "Max fluctuation in time between two signaling messages sent by an agent" , 0 , true , Double.MAX_VALUE , true);
079        private InputParameter signaling_averageDelay = new InputParameter ("signaling_averageDelay", 0.0 , "Average time between signaling message transmission by an agent and its reception by other or others" , 0 , true , Double.MAX_VALUE , true);
080        private InputParameter signaling_maxFluctuationInDelay = new InputParameter ("signaling_maxFluctuationInDelay", 0.0 , "Max fluctuation in time in the signaling delay, in absolute time values. The signaling delays are sampled from a uniform distribution within the given interval" , 0 , true , Double.MAX_VALUE , true);
081        private InputParameter signaling_signalingLossProbability = new InputParameter ("signaling_signalingLossProbability", 0.05 , "Probability that a signaling message transmitted is lost (not received by other or others involved agents)" , 0 , true , Double.MAX_VALUE , true);
082
083        private InputParameter routing_fixedPacketDurationAndSchedulingInterval = new InputParameter ("routing_fixedPacketDurationAndSchedulingInterval", 1.0 , "Fixed slot size (synchronized all links)" , 0 , false , Double.MAX_VALUE , true);
084        private InputParameter routing_numTrafficUnitsOfOnePacket = new InputParameter ("routing_numTrafficUnitsOfOnePacket", 1.0 , "One packet in one slot time is this traffic" , 0 , false , Double.MAX_VALUE , true);
085        private InputParameter routing_maxNumberOfN2PStatisticsIntervals = new InputParameter ("routing_maxNumberOfN2PStatisticsIntervals", (int) 50 , "Simulation length in number of scheduling intervals" , 1 , Integer.MAX_VALUE);
086        private InputParameter routing_statNumSchedSlotBetweenN2PRecomputing = new InputParameter ("routing_statNumSchedSlotBetweenN2PRecomputing", (int) 100 , "Number of scheduling slots between two N2P update and storing a stat trace" , 1 , Integer.MAX_VALUE);
087        private InputParameter routing_pressureDifferenceThreshold = new InputParameter ("routing_pressureDifferenceThreshold", 1.0 , "Pressure difference parameter (delta P in the book)" , 0 , false , Double.MAX_VALUE , true);
088
089        private InputParameter gradient_gammaStep = new InputParameter ("gradient_gammaStep", 0.01 , "Gamma step in the gradient algorithm" , 0 , false , Double.MAX_VALUE , true);
090
091        private InputParameter simulation_randomSeed = new InputParameter ("simulation_randomSeed", (long) 1 , "Seed of the random number generator");
092        private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "backpressureRoutingDual" , "Root of the file name to be used in the output files. If blank, no output");
093
094        private NetPlan currentNetPlan;
095
096        private Map<List<Link>,Route> stat_mapSeqLinks2RouteId;
097        private Map<Route,Integer> stat_mapRouteId2CarriedPacketsLastInterval;
098        private int [] stat_accumNumGeneratedPackets_d;
099        private int [] stat_accumNumReceivedPackets_d;
100        private int [][] stat_accumNumQeueusPackets_nd;
101        private double stat_totalOfferedTrafficConstant;
102        
103        private TimeTrace stat_traceOf_xp; // traces are taken from netPlan, when it is updated (stat_n2pRecomputingInterval)
104        private TimeTrace stat_traceOf_objFunction;
105        private TimeTrace stat_traceOf_ye;
106        private TimeTrace stat_traceOf_queueSizes;
107
108        private IntMatrix3D ctlMostUpdatedQndKnownByNodeN1_n1nd; // sparse for memory efficiency 
109        private int [][] ctlNumPacketsQueue_nd;
110        private Map<Pair<Node,Demand>,LinkedList<PacketInfo>> ctlQueue_nd; 
111
112        private NetPlan optNetPlan;
113        private double [][] optQueueSizes_nd;
114        
115        @Override
116        public String getDescription()
117        {
118                return "This module implements a distributed dual-gradient based algorithm for adapting the network routing to the one which minimizes the average number of hops, that results in a purely decentralized backpressure scheme.";
119        }
120
121        @Override
122        public List<Triple<String, String, String>> getParameters()
123        {
124                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
125                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
126        }
127
128        @Override
129        public void initialize(NetPlan currentNetPlan, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
130        {
131                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
132                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
133
134                if (currentNetPlan.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks");
135
136                this.rng = new Random (simulation_randomSeed.getLong());
137                this.currentNetPlan = currentNetPlan;
138
139                /* Update the indexes. The set of nodes, links and demands cannot change during the algorithm */
140                this.N = this.currentNetPlan.getNumberOfNodes();
141                this.E = this.currentNetPlan.getNumberOfLinks();
142                this.D = this.currentNetPlan.getNumberOfDemands();
143                if ((E == 0) || (D == 0)) throw new Net2PlanException ("The input design should have links and demands");
144                
145                /* Sets the initial routing, with all the traffic balanced equally in all the paths of the demand */
146                this.stat_mapSeqLinks2RouteId = new HashMap<List<Link>,Route> ();
147                this.stat_mapRouteId2CarriedPacketsLastInterval = new HashMap<Route,Integer> ();
148
149                this.currentNetPlan.removeAllUnicastRoutingInformation();
150                this.currentNetPlan.setRoutingType(RoutingType.SOURCE_ROUTING);
151                final Pair<NetPlan,double [][]> pOpt = computeOptimumSolution (false);
152                final Pair<NetPlan,double [][]> pOpt01 = computeOptimumSolution (true);
153                System.out.println("opt No 01" + DoubleFactory2D.dense.make(pOpt.getSecond()));
154                System.out.println("opt 01" + DoubleFactory2D.dense.make(pOpt01.getSecond()));
155
156//              if (1==1) throw new RuntimeException ("Bad");
157                
158                this.optNetPlan = pOpt.getFirst();
159                this.optQueueSizes_nd = pOpt.getSecond();
160                /* Add the routes in the optimum solution to the netPlan, with zero traffic. For having teh same IDs later */
161                this.currentNetPlan.copyFrom(optNetPlan);
162                for (Route r: currentNetPlan.getRoutes())
163                {
164                        r.setCarriedTraffic(0.0 , 0.0);
165                        this.stat_mapSeqLinks2RouteId.put(r.getSeqLinksRealPath(),r);
166                        this.stat_mapRouteId2CarriedPacketsLastInterval.put(r, 0);
167                }
168                this.ctlQueue_nd = new HashMap<Pair<Node,Demand>,LinkedList<PacketInfo>> ();
169                for (Node n: this.currentNetPlan.getNodes())
170                        for (Demand d: this.currentNetPlan.getDemands())
171                                { ctlQueue_nd.put(Pair.of(n, d), new LinkedList<PacketInfo> ());  }
172
173                this.ctlMostUpdatedQndKnownByNodeN1_n1nd = IntFactory3D.sparse.make(N, N, D); 
174                this.ctlNumPacketsQueue_nd = new int [N][D];
175                
176                this.stat_accumNumGeneratedPackets_d = new int [D];
177                this.stat_accumNumReceivedPackets_d = new int [D];
178                this.stat_accumNumQeueusPackets_nd= new int [N][D]; 
179                
180                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
181                this.scheduleEvent(new SimEvent (routing_statNumSchedSlotBetweenN2PRecomputing.getInt() * routing_fixedPacketDurationAndSchedulingInterval.getDouble(), SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_STATISTICTRACES , -1));
182                this.scheduleEvent(new SimEvent (routing_fixedPacketDurationAndSchedulingInterval.getDouble() , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_INITTIMESLOTFORSCHEDULINGPACKETS , -1));
183                for (Node n : currentNetPlan.getNodes())
184                {
185                        final double signalingTime = signaling_isSynchronous.getBoolean()? signaling_averageInterMessageTime.getDouble() : Math.max(0 , signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
186                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , n));
187                }
188
189                /* Intialize the traces */
190                this.stat_totalOfferedTrafficConstant = this.currentNetPlan.getVectorDemandOfferedTraffic().zSum();
191                this.stat_traceOf_xp = new TimeTrace ();
192                this.stat_traceOf_ye = new TimeTrace ();
193                this.stat_traceOf_queueSizes = new TimeTrace ();
194                this.stat_traceOf_objFunction = new TimeTrace (); 
195                this.stat_traceOf_xp.add(0.0 , netPlanRouteCarriedTrafficMap (this.currentNetPlan));
196                this.stat_traceOf_ye.add(0.0, this.currentNetPlan.getVectorLinkTotalCarriedTraffic());
197                this.stat_traceOf_queueSizes.add(0.0, copyOf(this.ctlNumPacketsQueue_nd));
198                this.stat_traceOf_objFunction.add(0.0, computeObjectiveFucntionFromNetPlan());
199                
200        }
201
202        @Override
203        public void processEvent(NetPlan currentNetPlan, SimEvent event)
204        {
205                final double time = event.getEventTime();
206                switch (event.getEventType())
207                {
208                case UPDATE_STATISTICTRACES:
209                {
210                        for (Route r : this.currentNetPlan.getRoutes())
211                        {
212                                final int numPacketsArrivedDestinationLastInterval = stat_mapRouteId2CarriedPacketsLastInterval.get(r);
213                                final double trafficArrivedDestinationLastInterval = numPacketsArrivedDestinationLastInterval * routing_numTrafficUnitsOfOnePacket.getDouble() / (routing_fixedPacketDurationAndSchedulingInterval.getDouble() * routing_statNumSchedSlotBetweenN2PRecomputing.getInt());
214                                r.setCarriedTraffic(trafficArrivedDestinationLastInterval , trafficArrivedDestinationLastInterval);
215
216                                /* Set zero in the counter of packets */
217                                stat_mapRouteId2CarriedPacketsLastInterval.put(r,0);
218                        }
219                        
220                        /* Update the traces */
221                        this.stat_traceOf_xp.add(time, netPlanRouteCarriedTrafficMap(this.currentNetPlan));
222                        this.stat_traceOf_ye.add(time, this.currentNetPlan.getVectorLinkTotalCarriedTraffic());
223                        this.stat_traceOf_objFunction.add(time, computeObjectiveFucntionFromNetPlan());
224                        final double scaleFactorAccumNumQueuePacketsToAverageQueuedTraffic = this.routing_numTrafficUnitsOfOnePacket.getDouble() / this.routing_statNumSchedSlotBetweenN2PRecomputing.getInt();
225                        /* We store the average queue sizes in traffic units */
226                        this.stat_traceOf_queueSizes.add(time, scaledCopyOf(this.stat_accumNumQeueusPackets_nd , scaleFactorAccumNumQueuePacketsToAverageQueuedTraffic));
227                        this.stat_accumNumQeueusPackets_nd = new int [N][D]; // reset the count
228                        
229                        this.scheduleEvent(new SimEvent (time + routing_fixedPacketDurationAndSchedulingInterval.getDouble() * routing_statNumSchedSlotBetweenN2PRecomputing.getInt() , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_STATISTICTRACES , -1));
230
231                        break;
232                }
233                
234                case SIGNALING_RECEIVEDMESSAGE: // A node receives from an out neighbor the q_nd for any demand
235                {
236                        final Triple<Node,Node,int []> signalInfo = (Triple<Node,Node,int []>) event.getEventObject();
237                        final Node nMe = signalInfo.getFirst();
238                        final int index_nIdMe = nMe.getIndex ();
239                        final Node nNeighbor = signalInfo.getSecond();
240                        final int index_nNeighbor = nNeighbor.getIndex ();
241                        final int [] receivedInfo_q_d = signalInfo.getThird();
242                        for (int index_d = 0 ; index_d < D ; index_d ++)
243                                ctlMostUpdatedQndKnownByNodeN1_n1nd.set(index_nIdMe, index_nNeighbor, index_d, receivedInfo_q_d[index_d]);
244                        break;
245
246                }
247                
248                case SIGNALING_WAKEUPTOSENDMESSAGE: // A node broadcasts signaling info to its 1 hop neighbors
249                {
250                        final Node nMe = (Node) event.getEventObject();
251                        final int index_nMe = nMe.getIndex ();
252                        
253                        log.println("t=" + time + ": ROUTING_NODEWAKEUPTO_TRANSMITSIGNAL node " + nMe);
254
255                        /* Create the info I will signal */
256                        int [] infoToSignal_q_d = Arrays.copyOf(ctlNumPacketsQueue_nd [index_nMe], D);
257
258                        /* Send the events of the signaling information messages to incoming neighbors */
259                        for (Node n : nMe.getInNeighbors())
260                        {
261                                if (rng.nextDouble() >= this.signaling_signalingLossProbability.getDouble()) // the signaling may be lost => lost to all nodes
262                                {
263                                        final double signalingReceptionTime = time + Math.max(0 , signaling_averageDelay.getDouble() + signaling_maxFluctuationInDelay.getDouble() * (rng.nextDouble() - 0.5));
264                                        this.scheduleEvent(new SimEvent (signalingReceptionTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_RECEIVEDMESSAGE , Triple.of(n , nMe , infoToSignal_q_d)));
265                                }
266                        }
267                        
268                        /* Re-schedule when to wake up again */ 
269                        final double signalingTime = signaling_isSynchronous.getBoolean()? time + signaling_averageInterMessageTime.getDouble() : Math.max(time , time + signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
270                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , nMe));
271                        break;
272                }
273
274                case UPDATE_INITTIMESLOTFORSCHEDULINGPACKETS: // a node updates its p_n, p_e values, using most updated info available
275                {
276                        log.println("t=" + time + ": UPDATE_INITTIMESLOTFORSCHEDULINGPACKETS");
277
278                        /* First, create traffic of the demands. A constant amount */
279                        for (Demand d : this.currentNetPlan.getDemands())
280                        {
281                                final double avgNumPackets_d = d.getOfferedTraffic() / routing_numTrafficUnitsOfOnePacket.getDouble();
282                                final int numPackets_d = Poisson.staticNextInt(avgNumPackets_d);
283                                final Node a_d = d.getIngressNode();
284                                LinkedList<PacketInfo> q = this.ctlQueue_nd.get(Pair.of(a_d, d));
285                                for (int cont = 0 ; cont < numPackets_d ; cont ++)
286                                        q.add (new PacketInfo (d , time));
287                                this.ctlNumPacketsQueue_nd [a_d.getIndex ()][d.getIndex ()] += numPackets_d;
288                                this.stat_accumNumGeneratedPackets_d [d.getIndex ()] += numPackets_d;
289                        }
290
291                        /* Forwarding decisions for each link. Each link can transmit as many packets as its capacity divided by trafic per packet factor (floor), during this scheduling */
292                        int [][] numPacketsToForward_ed = new int [E][D];
293                        int [][] numPacketsToForward_nd = new int [N][D];
294                        for (int index_e = 0 ; index_e  < E ; index_e ++)
295                        {
296                                final Link e = currentNetPlan.getLink (index_e);
297                                final int maxNumPacketsToTransmit = (int) Math.floor(e.getCapacity() / this.routing_numTrafficUnitsOfOnePacket.getDouble());
298                                final Node a_e = e.getOriginNode();
299                                final Node b_e = e.getDestinationNode();
300                                final int index_ae = a_e.getIndex ();
301                                final int index_be = b_e.getIndex ();
302                                for (int cont = 0 ; cont < maxNumPacketsToTransmit ; cont ++)
303                                {
304                                        /* Take the queue with higher */
305                                        double highestBackpressure = -1;
306                                        int indexDemandToTransmitPacket = -1;
307                                        for (int index_d = 0 ; index_d < D ; index_d ++)
308                                        {
309                                                final Demand d = currentNetPlan.getDemand (index_d);
310                                                final int thisQueueSize = ctlNumPacketsQueue_nd [index_ae][index_d];
311                                                if (thisQueueSize != ctlQueue_nd.get(Pair.of(a_e, d)).size ()) throw new RuntimeException ("Bad"); 
312                                                final int neighborQueueSize = ctlMostUpdatedQndKnownByNodeN1_n1nd.get(index_ae,index_be,index_d);
313                                                final double backpressure = gradient_gammaStep.getDouble() * (thisQueueSize - numPacketsToForward_nd[index_ae][index_d] - neighborQueueSize);
314                                                if (backpressure > highestBackpressure) { indexDemandToTransmitPacket = index_d; highestBackpressure = backpressure; } 
315                                        }
316                                        if (highestBackpressure <= routing_pressureDifferenceThreshold.getDouble()) break; // no more queues will be empty, not enough backpressures
317                                        
318                                        numPacketsToForward_ed[index_e][indexDemandToTransmitPacket] ++;
319                                        numPacketsToForward_nd[index_ae][indexDemandToTransmitPacket] ++;
320                                }
321                        }
322                        /* Apply the forwarding decisions */
323                        for (int index_e = 0 ; index_e  < E ; index_e ++)
324                        {
325                                final Link e = currentNetPlan.getLink (index_e);
326                                final Node a_e = e.getOriginNode();
327                                final Node b_e = e.getDestinationNode();
328                                final int index_ae = a_e.getIndex ();
329                                final int index_be = b_e.getIndex ();
330                                for (int index_d = 0 ; index_d < D ; index_d ++)
331                                {
332                                        final Demand d = currentNetPlan.getDemand (index_d);
333                                        final int numPackets_ed = numPacketsToForward_ed [index_e][index_d];
334
335                                        for (int cont = 0 ; cont < numPackets_ed ; cont ++)
336                                        {
337                                                PacketInfo p = this.ctlQueue_nd.get(Pair.of(a_e,d)).pop();
338                                                ctlNumPacketsQueue_nd [index_ae][index_d] --;
339                                                if (this.ctlQueue_nd.get(Pair.of(a_e,d)).size() != ctlNumPacketsQueue_nd [index_ae][index_d]) throw new RuntimeException ("Bad");
340                                                p.seqLinks.add(e); // update the seq of links
341                                                if (b_e == d.getEgressNode())
342                                                {
343                                                        /* Packet arrived to destination: update the statistics */
344                                                        Route route = stat_mapSeqLinks2RouteId.get(p.seqLinks);
345                                                        
346                                                        /* If the route does not exist, create it */
347                                                        if (route == null) 
348                                                        {
349                                                                route = this.currentNetPlan.addRoute(d, 0.0 , 0.0 , p.seqLinks, null);
350                                                                stat_mapSeqLinks2RouteId.put(p.seqLinks , route);
351                                                                stat_mapRouteId2CarriedPacketsLastInterval.put(route, 1);
352                                                        }
353                                                        else
354                                                        {
355                                                                stat_mapRouteId2CarriedPacketsLastInterval.put(route, stat_mapRouteId2CarriedPacketsLastInterval.get(route) + 1);
356                                                                this.stat_accumNumReceivedPackets_d [index_d] ++;
357                                                        }
358                                                }
359                                                else
360                                                {
361                                                        /* Not the end node => put it in the next queue */
362                                                        this.ctlQueue_nd.get(Pair.of(b_e,d)).addLast(p);
363                                                        ctlNumPacketsQueue_nd [index_be][index_d] ++;
364                                                }
365                                        }
366                                }
367                                
368                        }
369                        
370                        /* Update the statistics of accum number of packets in a queue */
371                        for (int index_n = 0 ; index_n < N ; index_n ++)
372                                for (int index_d = 0 ; index_d < D ; index_d ++)
373                                        stat_accumNumQeueusPackets_nd [index_n][index_d] += this.ctlNumPacketsQueue_nd [index_n][index_d];
374                        
375
376                        this.scheduleEvent(new SimEvent (time + routing_fixedPacketDurationAndSchedulingInterval.getDouble() , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_INITTIMESLOTFORSCHEDULINGPACKETS , -1));
377        
378                        if (time > this.routing_maxNumberOfN2PStatisticsIntervals.getInt() * routing_statNumSchedSlotBetweenN2PRecomputing.getInt() * this.routing_fixedPacketDurationAndSchedulingInterval.getDouble()) { this.endSimulation (); }
379                        break;
380                }
381                        
382
383                default: throw new RuntimeException ("Unexpected received event");
384                }
385                
386                
387        }
388
389        public String finish (StringBuilder st , double simTime)
390        {
391                System.out.println("stat_accumNumGeneratedPackets_d: " + Arrays.toString(stat_accumNumGeneratedPackets_d));
392                System.out.println("stat_accumNumReceivedPackets_d: " + Arrays.toString(stat_accumNumReceivedPackets_d));
393                System.out.println("transmittedNotReceived: " + Arrays.toString(IntUtils.substract(stat_accumNumGeneratedPackets_d, stat_accumNumReceivedPackets_d)));
394                double [] avTrafficGenerated_d = new double [D];
395                for (int index_d = 0; index_d < D ; index_d ++)
396                        avTrafficGenerated_d[index_d] = stat_accumNumGeneratedPackets_d [index_d] * routing_numTrafficUnitsOfOnePacket.getDouble() / simTime;
397
398//              System.out.println("avTrafficGenerated_d: " + Arrays.toString(avTrafficGenerated_d));
399//              System.out.println("avTrafficGenerated_d - offeredFromN2P: " + Arrays.toString(DoubleUtils.substract(avTrafficGenerated_d, this.currentNetPlan.getVectorDemandOfferedTraffic().toArray())));
400
401                /* If no output file, return */
402                if (simulation_outFileNameRoot.getString().equals("")) return null;
403                
404                /* Compute optimum solution and cost */
405                double optCost = optNetPlan.getVectorLinkTotalCarriedTraffic().zSum() / stat_totalOfferedTrafficConstant;
406                
407                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_objFunc.txt"), optCost);
408                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_xp.txt"), optNetPlan.getVectorRouteCarriedTraffic());
409                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_ye.txt"), optNetPlan.getVectorLinkTotalCarriedTraffic());
410                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_qnd.txt"), optQueueSizes_nd);
411                this.stat_traceOf_queueSizes.printToFile(new File (simulation_outFileNameRoot.getString() + "_qnd.txt"));
412                this.stat_traceOf_objFunction.printToFile(new File (simulation_outFileNameRoot.getString() + "_objFunc.txt"));
413                this.stat_traceOf_ye.printToFile(new File (simulation_outFileNameRoot.getString() + "_ye.txt"));
414                for (Pair<Double,Object> sample : this.stat_traceOf_xp.getList())
415                {
416                        Map<Long,Double> x_p = (Map<Long,Double>) sample.getSecond ();
417                        for (long r: this.currentNetPlan.getRouteIds())
418                                if (!x_p.containsKey(r)) x_p.put(r , 0.0); 
419                }
420                this.stat_traceOf_xp.printToFile(new File (simulation_outFileNameRoot.getString() + "_xp.txt"));
421                return null;
422        }
423        
424        private double computeObjectiveFucntionFromNetPlan ()
425        {
426                return this.currentNetPlan.getVectorLinkTotalCarriedTraffic().zSum () / this.stat_totalOfferedTrafficConstant;
427        }
428
429        private Pair<NetPlan,double [][]> computeOptimumSolution (boolean xdeVariablesAsFractionsOfTraffic)
430        {
431                /* Modify the map so that it is the pojection where all elements sum h_d, and are non-negative */
432                final int D = this.currentNetPlan.getNumberOfDemands();
433                final int E = this.currentNetPlan.getNumberOfLinks();
434                final DoubleMatrix1D h_d = this.currentNetPlan.getVectorDemandOfferedTraffic();
435                final DoubleMatrix1D u_e = this.currentNetPlan.getVectorLinkCapacity();
436                                
437                OptimizationProblem op = new OptimizationProblem();
438                if (xdeVariablesAsFractionsOfTraffic)
439                        op.addDecisionVariable("x_de", false , new int[] { D , E }, 0 , 1);
440                else
441                        op.addDecisionVariable("x_de", false , new int[] { D , E }, 0 , u_e.getMaxLocation() [0]);
442                op.setInputParameter("h_d", h_d , "row");
443                op.setInputParameter("u_e", u_e , "row");
444                op.setInputParameter("deltaP", routing_pressureDifferenceThreshold.getDouble());
445                
446                /* Set some input parameters */
447                if (xdeVariablesAsFractionsOfTraffic)
448                        op.setObjectiveFunction("minimize", "deltaP * sum (h_d * x_de)");
449                else
450                        op.setObjectiveFunction("minimize", "deltaP * sum (x_de)");
451
452                /* VECTORIAL FORM OF THE CONSTRAINTS  */
453                op.setInputParameter("A_ne", this.currentNetPlan.getMatrixNodeLinkIncidence());
454                op.setInputParameter("A_nd", this.currentNetPlan.getMatrixNodeDemandIncidence());
455                if (xdeVariablesAsFractionsOfTraffic)
456                {
457                        op.addConstraint("A_ne * (x_de') >= A_nd" , "flowConservationConstraints"); // the flow-conservation constraints (NxD constraints)
458                        op.addConstraint("h_d * x_de <= u_e"); // the capacity constraints (E constraints)
459                }
460                else
461                {
462                        op.addConstraint("A_ne * (x_de') >= A_nd * diag(h_d)" , "flowConservationConstraints"); // the flow-conservation constraints (NxD constraints)
463                        op.addConstraint("sum(x_de,1) <= u_e"); // the capacity constraints (E constraints)
464                }
465                
466                /* Call the solver to solve the problem */
467                op.solve("cplex");
468
469                /* If an optimal solution was not found, quit */
470                if (!op.solutionIsOptimal ()) throw new RuntimeException ("An optimal solution was not found");
471        
472                /* Retrieve the optimum solutions. Convert the bps into Erlangs */
473                final DoubleMatrix2D x_de_array = op.getPrimalSolution("x_de").view2D ();
474                final double [][] q_nd_array = (double [][]) op.getMultipliersOfConstraint("flowConservationConstraints").toArray();
475                
476                //System.out.println("x_de:" + op.getPrimalSolution("x_de"));
477                
478                /* Convert the x_de variables into a set of routes for each demand  */
479                NetPlan np = this.currentNetPlan.copy();
480                np.removeAllUnicastRoutingInformation();
481                np.setRoutingFromDemandLinkCarriedTraffic(x_de_array , xdeVariablesAsFractionsOfTraffic , false);
482                
483                /* Check solution: all traffic is carried, no link oversubscribed */
484                for (Demand d : np.getDemands()) if (d.getBlockedTraffic() > 1E-3) throw new RuntimeException ("d: " + d + ", hd: " +  d.getOfferedTraffic() + ", carried_d: " + d.getCarriedTraffic() + "... Bad");
485                for (Link e : np.getLinks()) if (e.getCarriedTrafficIncludingProtectionSegments() > e.getCapacity() + 1E-3) throw new RuntimeException ("Bad");
486
487                /* Scale q_nd multipliers by 1/gamma factor => they become now the optimum queue sizes in the algorithm */
488                double [][] optQueueSizes = new double [N][D];
489                for (int i = 0 ; i < N ; i ++)
490                        for (int j = 0 ; j < D ; j ++)
491                                optQueueSizes [i][j] = q_nd_array [i][j] / this.gradient_gammaStep.getDouble();
492
493                return Pair.of(np,optQueueSizes);
494        }
495
496        private Map<Long,Double> netPlanRouteCarriedTrafficMap (NetPlan np) { Map<Long,Double> res = new HashMap<Long,Double> (); for (Route r : np.getRoutes ()) res.put (r.getId () , r.getCarriedTraffic()); return res; }
497        private class PacketInfo
498        {
499                final Demand d;
500                final double packetCreationTime; 
501                LinkedList<Link> seqLinks;
502                PacketInfo (Demand d , double creationTime) { this.d = d; this.packetCreationTime = creationTime; this.seqLinks = new LinkedList<Link> (); }
503        }
504        private int [][] copyOf (int [][] a) { int [][] res = new int [a.length][]; for (int cont = 0 ; cont < a.length ; cont ++) res [cont] = Arrays.copyOf(a[cont], a[cont].length); return res; }
505        private double [][] copyOf (double [][] a) { double [][] res = new double [a.length][]; for (int cont = 0 ; cont < a.length ; cont ++) res [cont] = Arrays.copyOf(a[cont], a[cont].length); return res; }
506        private double [][] scaledCopyOf (int [][] a , double mult) { double [][] res = new double [a.length][a[0].length]; for (int cont = 0 ; cont < a.length ; cont ++) for (int cont2 = 0; cont2 < a[cont].length ; cont2 ++) res [cont][cont2] = a[cont][cont2] * mult; return res; }
507}