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
013
014import java.io.File;
015import java.io.FileOutputStream;
016import java.io.PrintStream;
017import java.util.HashMap;
018import java.util.List;
019import java.util.Map;
020import java.util.Map.Entry;
021import java.util.Random;
022import java.util.Set;
023
024import cern.colt.matrix.tdouble.DoubleFactory1D;
025import cern.colt.matrix.tdouble.DoubleFactory2D;
026import cern.colt.matrix.tdouble.DoubleMatrix1D;
027import cern.colt.matrix.tdouble.DoubleMatrix2D;
028
029import com.net2plan.examples.ocnbook.offline.Offline_ca_wirelessPersistenceProbability;
030import com.net2plan.interfaces.networkDesign.Configuration;
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.simulation.IEventProcessor;
036import com.net2plan.interfaces.simulation.SimEvent;
037import com.net2plan.libraries.NetworkPerformanceMetrics;
038import com.net2plan.libraries.WirelessUtils;
039import com.net2plan.utils.GradientProjectionUtils;
040import com.net2plan.utils.InputParameter;
041import com.net2plan.utils.TimeTrace;
042import com.net2plan.utils.Triple;
043
044/** 
045 * This module implements a distributed primal-gradient based algorithm for adjusting the link persistence probabilities in a wireless network with a ALOHA-type random-access based MAC, to maximize the network utility enforcing a fair allocation of the resources.
046 *
047 * Ths event processor is adapted to permit observing the algorithm performances under user-defined conditions, 
048 * including asynchronous distributed executions, where signaling can be affected by losses and/or delays, and/or measurement errors. 
049 * The time evolution of different metrics can be stored in output files, for later processing. 
050 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec9_5_persitenceProbabilityAdjustmentPrimal.m">{@code fig_sec9_5_persitenceProbabilityAdjustmentPrimal.m}</a> MATLAB file used for generating the graph/s of the case study in the 
051 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm.
052 * 
053 * To simulate a network with this module, use the {@code Online_evGen_doNothing} generator.
054 * 
055 * @net2plan.keywords Random-access MAC, Wireless, Distributed algorithm, Primal gradient algorithm, Capacity assignment (CA)
056 * @net2plan.ocnbooksections Section 9.5
057 * @net2plan.inputParameters 
058 * @author Pablo Pavon-Marino
059 */
060public class Online_evProc_persistenceProbAdjustmentPrimal extends IEventProcessor
061{
062        private static PrintStream getNulFile () { try { return new PrintStream (new FileOutputStream ("NUL") , false); } catch (Exception e) {e.printStackTrace(); throw new RuntimeException ("Not NUL file"); }   } 
063        private PrintStream log = getNulFile (); //System.err;//new PrintStream (new FileOutputStream ("NUL") , true);
064
065        private Random rng;
066        
067        private static final int SIGNALING_WAKEUPTOSENDMESSAGE = 300;
068        private static final int SIGNALING_RECEIVEDMESSAGE = 301;
069        private static final int UPDATE_WAKEUPTOUPDATE = 302;
070
071        private InputParameter signaling_isSynchronous = new InputParameter ("signaling_isSynchronous", false , "true if all the distributed agents involved wake up synchronously to send the signaling messages");
072        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);
073        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);
074        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);
075        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);
076        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);
077        private InputParameter update_isSynchronous = new InputParameter ("update_isSynchronous", false , "true if all the distributed agents involved wake up synchronousely to update its state");
078        private InputParameter update_averageInterUpdateTime = new InputParameter ("update_averageInterUpdateTime", 1.0 , "Average time between two updates of an agent" , 0 , false , Double.MAX_VALUE , true);
079        private InputParameter update_maxFluctuationInterUpdateTime = new InputParameter ("update_maxFluctuationInterUpdateTime", 0.5 , "Max fluctuation in time in the update interval of an agent, in absolute time values. The update intervals are sampled from a uniform distribution within the given interval" , 0 , true , Double.MAX_VALUE , true);
080
081        private InputParameter gradient_maxGradientAbsoluteNoise = new InputParameter ("gradient_maxGradientAbsoluteNoise", 0.0 , "Max value of the added noise to the gradient coordinate in absolute values" , 0 , true , Double.MAX_VALUE , true);
082        private InputParameter gradient_gammaStep = new InputParameter ("gradient_gammaStep", 0.00001 , "Gamma step in the gradient algorithm" , 0 , false , Double.MAX_VALUE , true);
083        private InputParameter gradient_maxGradientCoordinateChange = new InputParameter ("gradient_maxGradientCoordinateChange", 0.1 , "Maximum change in an iteration of a gradient coordinate" , 0 , false , Double.MAX_VALUE , true);
084        private InputParameter simulation_randomSeed = new InputParameter ("simulation_randomSeed", (long) 1 , "Seed of the random number generator");
085        private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "persistenceProbAdjustmentPrimal" , "Root of the file name to be used in the output files. If blank, no output");
086        private InputParameter simulation_maxNumberOfUpdateIntervals = new InputParameter ("simulation_maxNumberOfUpdateIntervals", 100.0 , "Maximum number of update intervals in average per agent" , 0 , false , Double.MAX_VALUE , true);
087
088        private InputParameter control_maxUeRelativeNoise = new InputParameter ("control_maxUeRelativeNoise", 0 , "Maximum relative fluctuation in link capacity estimation, that reflects into gradient noise" , 0 , true , Double.MAX_VALUE , true);
089        private InputParameter control_simultaneousTxAndRxPossible = new InputParameter ("control_simultaneousTxAndRxPossible", false , "If false, a node cannot transmit and receive simultaneously. If true, it can. This affects the interference map calculation");
090        private InputParameter control_linkNominalCapacity = new InputParameter ("control_linkNominalCapacity", 1.0 , "Nominal rate of the links. If non positive, nominal rates are rates of the initial design");
091        private InputParameter control_minLinkPersistenceProb = new InputParameter ("control_minLinkPersistenceProb", 0.03 , "Minimum persistence probability of a link" , 0 , true , 1 , true);
092        private InputParameter control_maxNodePersistenceProb = new InputParameter ("control_maxNodePersistenceProb", 0.99 , "Maximum persistence probability of a node (summing the persistence probabilities of its outoingl links)" , 0 , true , 1 , true);
093        private InputParameter control_fairnessFactor = new InputParameter ("control_fairnessFactor", 2.0 , "Fairness factor in utility function of link capacities (>= 1 to be a convex program)" , 0 , true , Double.MAX_VALUE , true);
094
095        private NetPlan currentNetPlan , copyInitialNetPlan;
096        private int E , N;
097                
098        private DoubleMatrix1D control_linkNominalCapacities;
099        private DoubleMatrix1D control_p_e;
100        private DoubleMatrix1D control_q_n;
101        private DoubleMatrix2D control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2; // each node is assigned a map with known p_n values of one hop neighbors and two hop neighbors
102        private DoubleMatrix2D control_mostUpdatedUeValueNodeKnowByNode_ne; // each node is assigned a map with known s_e values of one hop neighbors and two hop neighbors
103
104        private Map<Link,Set<Node>> control_nodesInterfTo_e;
105        private Map<Node,Set<Link>> control_linksInterfFrom_n;
106        private Map<Node,int []> control_outLinksIndexes;
107        
108        private TimeTrace stat_traceOf_p_e;
109        private TimeTrace stat_traceOf_u_e;
110        private TimeTrace stat_traceOf_objFunction;
111
112        
113        @Override
114        public String getDescription()
115        {
116                return "This module implements a distributed primal-gradient based algorithm for adjusting the link persistence probabilities in a wireless network with a ALOHA-type random-access based MAC, to maximize the network utility enforcing a fair allocation of the resources.";
117        }
118
119        @Override
120        public List<Triple<String, String, String>> getParameters()
121        {
122                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
123                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
124        }
125
126        @Override
127        public void initialize(NetPlan currentNetPlan, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
128        {
129                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
130                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
131
132                this.currentNetPlan = currentNetPlan;
133                this.copyInitialNetPlan = currentNetPlan.copy ();
134                if (currentNetPlan.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks");
135                this.N = currentNetPlan.getNumberOfNodes ();
136                this.E = currentNetPlan.getNumberOfLinks ();
137
138                if (control_linkNominalCapacity.getDouble() > 0) currentNetPlan.setVectorLinkCapacity(DoubleFactory1D.dense.make (E , control_linkNominalCapacity.getDouble()));
139                control_linkNominalCapacities = currentNetPlan.getVectorLinkCapacity();
140                
141                this.rng = new Random (simulation_randomSeed.getLong());
142
143                /* Check all nodes have output links */
144                for (Node n : currentNetPlan.getNodes()) if (n.getOutgoingLinks().isEmpty()) throw new Net2PlanException ("All nodes should have output links");
145                for (Node n : currentNetPlan.getNodes()) if (control_minLinkPersistenceProb.getDouble() * n.getOutgoingLinks().size () > control_maxNodePersistenceProb.getDouble()) throw new Net2PlanException("Minimum persistence per link is too high or maximum persistence per node too small: the problem has no feasible solutions");
146
147                Triple<DoubleMatrix2D, Map<Link,Set<Node>> , Map<Node,Set<Link>>> interfMap = WirelessUtils.computeInterferenceMap(currentNetPlan.getNodes () , currentNetPlan.getLinks () , control_simultaneousTxAndRxPossible.getBoolean());
148                this.control_nodesInterfTo_e = interfMap.getSecond ();
149                this.control_linksInterfFrom_n = interfMap.getThird();
150                
151                this.control_outLinksIndexes = new HashMap<Node,int []> ();
152                for (Node n : currentNetPlan.getNodes ()) 
153                { 
154                        int [] indexes = new int [n.getOutgoingLinks().size ()]; 
155                        int counter = 0; for (Link e : n.getOutgoingLinks()) indexes [counter ++] = e.getIndex (); 
156                        control_outLinksIndexes.put (n , indexes);
157                }
158                
159                
160                /* Initialize the persistence probabilities: p_e to the minimum */
161                this.control_q_n = DoubleFactory1D.dense.make (N);
162                this.control_p_e = DoubleFactory1D.dense.make (E , control_minLinkPersistenceProb.getDouble());
163                for (Node n : currentNetPlan.getNodes ()) 
164                        this.control_q_n.set(n.getIndex () , control_minLinkPersistenceProb.getDouble() *  n.getOutgoingLinks().size ());
165                control_mostUpdatedUeValueNodeKnowByNode_ne = DoubleFactory2D.dense.make (N,E);
166                control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2 = DoubleFactory2D.dense.make (N,N);
167
168                /* Initialize the most updated signaled values known of q_n: one and two hop neighbors */
169                for (Node n1 : currentNetPlan.getNodes())
170                { 
171                        for (Node nNeighborOneHop : n1.getOutNeighbors())
172                        {
173                                control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2.set (n1.getIndex () , nNeighborOneHop.getIndex () , control_q_n.get (nNeighborOneHop.getIndex ()));
174                                for (Node nNeighborTwoHops : nNeighborOneHop.getOutNeighbors())
175                                        if (nNeighborTwoHops != n1) 
176                                                control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2.set (n1.getIndex () , nNeighborTwoHops.getIndex () , control_q_n.get (nNeighborTwoHops.getIndex ()));
177                        }
178                }
179                
180                /* Set link capacities to its actual values */
181                for (Link e : this.currentNetPlan.getLinks ()) 
182                {
183                        /* compute link capacity */
184                        double u_e = this.control_linkNominalCapacities.get(e.getIndex ()) * this.control_p_e.get(e.getIndex ());
185                        for (Node n : this.control_nodesInterfTo_e.get(e))
186                                u_e *= Math.max(0 , 1 - this.control_q_n.get(n.getIndex ()));
187                        e.setAttribute("p_e" , "" + control_p_e.get(e.getIndex ())); 
188                        e.setCapacity(u_e); 
189                }               
190
191                /* Set the most updated values known in each node, from the true u_e values */
192                DoubleMatrix1D noisyUeEstimation_e = DoubleFactory1D.dense.make (E);
193                for (Link e : currentNetPlan.getLinks())
194                        noisyUeEstimation_e.set(e.getIndex (), e.getCapacity() * ( 1 + control_maxUeRelativeNoise.getDouble()*2*(rng.nextDouble()-0.5)));
195                for (Node n : currentNetPlan.getNodes())
196                { 
197                        for (Link e : this.control_linksInterfFrom_n.get(n)) control_mostUpdatedUeValueNodeKnowByNode_ne.set (n.getIndex () , e.getIndex () , noisyUeEstimation_e.get(e.getIndex ()));
198                        for (Link e : n.getOutgoingLinks()) control_mostUpdatedUeValueNodeKnowByNode_ne.set (n.getIndex () , e.getIndex () , noisyUeEstimation_e.get(e.getIndex ()));
199                        for (Link e : n.getIncomingLinks()) control_mostUpdatedUeValueNodeKnowByNode_ne.set (n.getIndex () , e.getIndex () , noisyUeEstimation_e.get(e.getIndex ()));
200                }
201                
202                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
203                for (Node n : currentNetPlan.getNodes())
204                {
205                        final double signalingTime = (signaling_isSynchronous.getBoolean())? signaling_averageInterMessageTime.getDouble() : Math.max(0 , signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
206                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , n));
207                        final double updateTime = (update_isSynchronous.getBoolean())? update_averageInterUpdateTime.getDouble() : Math.max(0 , update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
208                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE , n));
209                }
210
211                /* Intialize the traces */
212                this.stat_traceOf_p_e = new TimeTrace ();
213                this.stat_traceOf_u_e = new TimeTrace ();
214                this.stat_traceOf_objFunction = new TimeTrace ();
215                this.stat_traceOf_p_e.add(0.0, control_p_e.copy());
216                this.stat_traceOf_u_e.add(0.0, currentNetPlan.getVectorLinkCapacity());
217                this.stat_traceOf_objFunction.add(0.0 , NetworkPerformanceMetrics.alphaUtility(currentNetPlan.getVectorLinkCapacity() , control_fairnessFactor.getDouble()));
218        }
219
220        @Override
221        public void processEvent(NetPlan currentNetPlan, SimEvent event)
222        {
223                final double t = event.getEventTime();
224                switch (event.getEventType())
225                {
226                case SIGNALING_RECEIVEDMESSAGE: // A node broadcasts signaling info to its 1 hop neighbors
227                {
228                        final Triple<Node,Map<Node,Double>,Map<Link,Double>> signalInfo = (Triple<Node,Map<Node,Double>,Map<Link,Double>>) event.getEventObject();
229                        final Node nMe = signalInfo.getFirst();
230                        final Map<Node,Double> receivedInfo_qn = signalInfo.getSecond();
231                        final Map<Link,Double> receivedInfo_ue = signalInfo.getThird();
232                        /* Take the p_n values signaled except my own p_n (I prefer my own estimation of q_n, than something estimated by others) */
233                        for (Entry<Node,Double> nodeInfo : receivedInfo_qn.entrySet())
234                                if (nodeInfo.getKey() != nMe) control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2.set (nMe.getIndex () , nodeInfo.getKey().getIndex () , nodeInfo.getValue());
235                        /* Take the u_e values signaled */
236                        for (Entry<Link,Double> linkInfo : receivedInfo_ue.entrySet())
237                                control_mostUpdatedUeValueNodeKnowByNode_ne.set (nMe.getIndex () , linkInfo.getKey().getIndex () , linkInfo.getValue());
238                        break;
239                }
240                
241                case SIGNALING_WAKEUPTOSENDMESSAGE: // A node broadcasts signaling info to its 1 hop neighbors
242                {
243                        final Node nMe = (Node) event.getEventObject();
244                        final DoubleMatrix1D infoIKnow_qn = control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2.viewRow (nMe.getIndex ());
245                        infoIKnow_qn.set (nMe.getIndex () , control_q_n.get (nMe.getIndex ()));
246
247                        /* Create the info I will signal */
248                        Map<Node,Double> infoToSignal_qn = new HashMap<Node,Double> ();
249                        infoToSignal_qn.put(nMe, this.control_q_n.get (nMe.getIndex ()));
250                        for (Node inNeighbor : nMe.getInNeighbors()) 
251                                infoToSignal_qn.put(inNeighbor, infoIKnow_qn.get(inNeighbor.getIndex ()));
252                        Map<Link,Double> infoToSignal_ue = new HashMap<Link,Double> ();
253                        for (Link inLink : nMe.getIncomingLinks()) 
254                                infoToSignal_ue.put(inLink, inLink.getCapacity() * ( 1 + 2*this.control_maxUeRelativeNoise.getDouble()*(rng.nextDouble() - 0.5)  ) );
255
256                        /* Send the events of the signaling information messages to the neighbors */
257                        if (rng.nextDouble() >= this.signaling_signalingLossProbability.getDouble()) // the signaling may be lost => lost to all nodes
258                                for (Node outNeighbor : nMe.getOutNeighbors())
259                                {
260                                        final double signalingReceptionTime = t + Math.max(0 , signaling_averageDelay.getDouble() + signaling_maxFluctuationInDelay.getDouble() * (rng.nextDouble() - 0.5));
261                                        this.scheduleEvent(new SimEvent (signalingReceptionTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_RECEIVEDMESSAGE , Triple.of(outNeighbor, infoToSignal_qn , infoToSignal_ue)));
262                                }
263                        
264                        /* Re-schedule when to wake up again */
265                        final double signalingTime = signaling_isSynchronous.getBoolean()? t + signaling_averageInterMessageTime.getDouble() : Math.max(t , t + signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
266                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , nMe));
267                        break;
268                }
269
270                case UPDATE_WAKEUPTOUPDATE: // a node updates its p_n, p_e values, using most updated info available
271                {
272                        final Node nMe = (Node) event.getEventObject(); 
273                        
274                        final DoubleMatrix1D infoIKnow_qn = control_mostUpdatedQn2ValueNodeKnowByNode1_n1n2.viewRow (nMe.getIndex ());
275                        infoIKnow_qn.set (nMe.getIndex () , control_q_n.get (nMe.getIndex ()));
276                        DoubleMatrix1D infoIKnow_ue = control_mostUpdatedUeValueNodeKnowByNode_ne.viewRow (nMe.getIndex ());
277                        DoubleMatrix1D newMacValues_e = DoubleFactory1D.dense.make (E);
278                        for (int eIndex : control_outLinksIndexes.get (nMe))
279                        {
280                                final Link e = currentNetPlan.getLink (eIndex);
281                                final double current_pe = this.control_p_e.get(eIndex);
282                                final double gradientThisLink = computeGradient (e , this.control_fairnessFactor.getDouble() , current_pe , infoIKnow_qn , infoIKnow_ue) + 2*gradient_maxGradientAbsoluteNoise.getDouble()*(rng.nextDouble()-0.5);
283                                newMacValues_e.set(eIndex, current_pe + gradient_gammaStep.getDouble() * gradientThisLink);
284                        }
285                        GradientProjectionUtils.euclideanProjection_sumInequality(newMacValues_e , control_outLinksIndexes.get (nMe) , control_minLinkPersistenceProb.getDouble() , control_maxNodePersistenceProb.getDouble());
286                        
287                        if (gradient_maxGradientCoordinateChange.getDouble() > 0)
288                                GradientProjectionUtils.scaleDown_maxAbsoluteCoordinateChange (control_p_e , newMacValues_e , control_outLinksIndexes.get (nMe) ,  gradient_maxGradientCoordinateChange.getDouble());
289                        for (int eIndex : control_outLinksIndexes.get (nMe)) this.control_p_e.set (eIndex , newMacValues_e.get (eIndex));
290                        
291                        /* update p_n value applied by me */
292                        double thisNode_pn = 0; for (Link e : nMe.getOutgoingLinks()) thisNode_pn += this.control_p_e.get(e.getIndex ());
293                        if ((thisNode_pn <= 0) || (thisNode_pn > 1 + 1E-3)) throw new RuntimeException ("Bad");
294                        this.control_q_n.set(nMe.getIndex (), thisNode_pn);
295                        
296                        /* Send event next recomputing time */
297                        final double updateTime = update_isSynchronous.getBoolean()? t + update_averageInterUpdateTime.getDouble() : Math.max(t , t + update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
298                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE,  nMe));
299
300                        /* Set link capacities to its actual values, update n2p object */
301                        for (Link e : currentNetPlan.getLinks ())
302                        { 
303                                /* compute link capacity */
304                                double u_e = this.control_linkNominalCapacities.get(e.getIndex ()) * this.control_p_e.get(e.getIndex ());
305                                for (Node n : this.control_nodesInterfTo_e.get(e))
306                                        u_e *= Math.max(0 , 1 - this.control_q_n.get(n.getIndex ()));
307                                e.setAttribute("p_e" , "" + control_p_e.get(e.getIndex ())); 
308                                e.setCapacity(u_e); 
309                        }               
310                        for (Node n : currentNetPlan.getNodes ()) { n.setAttribute("q_n" , "" + control_q_n.get(n.getIndex ())); }              
311
312                        checkCapacitiesNetPlan ();
313                        
314                        this.stat_traceOf_p_e.add(t, control_p_e.copy ());
315                        this.stat_traceOf_u_e.add(t, this.currentNetPlan.getVectorLinkCapacity());
316                        this.stat_traceOf_objFunction.add(t , NetworkPerformanceMetrics.alphaUtility(currentNetPlan.getVectorLinkCapacity() , control_fairnessFactor.getDouble()));
317
318                        if (t > this.simulation_maxNumberOfUpdateIntervals.getDouble() * this.update_averageInterUpdateTime.getDouble()) { this.endSimulation (); }
319                        
320                        break;
321                }
322                        
323
324                default: throw new RuntimeException ("Unexpected received event");
325                }
326                
327                
328        }
329
330        public String finish (StringBuilder st , double simTime)
331        {
332                if (simulation_outFileNameRoot.getString().equals("")) return null;
333                stat_traceOf_u_e.printToFile(new File (simulation_outFileNameRoot.getString() + "_ue.txt"));
334                stat_traceOf_objFunction.printToFile(new File (simulation_outFileNameRoot.getString() + "_objFunc.txt"));
335                stat_traceOf_p_e.printToFile(new File (simulation_outFileNameRoot.getString() + "_pe.txt"));
336                /* compute optimum solution */
337                Map<String,String> param = new HashMap<String,String> ();
338                param.put("alphaFairnessFactor", "" + this.control_fairnessFactor.getDouble());
339                param.put("minLinkPersistenceProb", "" + this.control_minLinkPersistenceProb.getDouble());
340                param.put("maxNodePersistenceProb", "" + this.control_maxNodePersistenceProb.getDouble());
341                param.put("linkNominalCapacity", "" + this.control_linkNominalCapacity.getDouble());
342                param.put("simultaneousTxAndRxPossible", "" + this.control_simultaneousTxAndRxPossible.getBoolean());
343                param.put("solverName", "ipopt");
344                param.put("solverLibraryName", "");
345                param.put("maxSolverTimeInSeconds", "-1");
346//              private InputParameter solverName = new InputParameter ("solverName", "#select# ipopt", "The solver name to be used by JOM. ");
347//              private InputParameter solverLibraryName = new InputParameter ("solverLibraryName", Configuration.getOption("glpkSolverLibraryName") , "The solver library full or relative path, to be used by JOM. Leave blank to use JOM default.");
348//              private InputParameter maxSolverTimeInSeconds = new InputParameter ("maxSolverTimeInSeconds", (double) -1 , "Maximum time granted to the solver to solve the problem. If this time expires, the solver returns the best solution found so far (if a feasible solution is found)");
349                new Offline_ca_wirelessPersistenceProbability ().executeAlgorithm(copyInitialNetPlan , param , null);
350                final double optimumNetUtilityJOM = NetworkPerformanceMetrics.alphaUtility(copyInitialNetPlan.getVectorLinkCapacity() , control_fairnessFactor.getDouble());
351                DoubleMatrix1D optimum_ue = copyInitialNetPlan.getVectorLinkCapacity();
352                DoubleMatrix1D optimum_pe = DoubleFactory1D.dense.make (E);
353                for (Link e : copyInitialNetPlan.getLinks ()) optimum_pe.set (e.getIndex () , Double.parseDouble (e.getAttribute("p_e")));
354                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_ue.txt"), optimum_ue);
355                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_pe.txt"), optimum_pe);
356                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_objFunc.txt"), optimumNetUtilityJOM);
357                return null;
358        }
359        
360        private double computeGradient (Link thisLink , double fairnessFactor , double thisLink_pe , DoubleMatrix1D infoIKnow_p_n , DoubleMatrix1D infoIKnow_ue)
361        {
362                final Node a_e = thisLink.getOriginNode();
363                if (thisLink_pe == 0) throw new RuntimeException ("Link " + thisLink + " has zero p_e");
364                double gradient = Math.pow(infoIKnow_ue.get(thisLink.getIndex ()) , 1 - fairnessFactor) / thisLink_pe; 
365                for (Link e : this.control_linksInterfFrom_n.get(a_e))
366                {
367                        if (e.getOriginNode() == thisLink.getOriginNode()) throw new RuntimeException ("Bad");
368                        gradient -= Math.pow(infoIKnow_ue.get(e.getIndex ()),1-fairnessFactor) / (1 - infoIKnow_p_n.get(e.getOriginNode().getIndex ()));
369                }
370                if (Double.isNaN(gradient)) throw new RuntimeException ("Link " + thisLink + " gradient is Nan");
371                if (Double.isInfinite(gradient)) throw new RuntimeException ("Link " + thisLink + " gradient is inifinte");
372                return gradient;
373        }
374
375        private void checkCapacitiesNetPlan ()
376        {
377                double [] p_e = new double [E]; 
378                for (Link e : currentNetPlan.getLinks ()) { p_e [e.getIndex ()] = Double.parseDouble (e.getAttribute("p_e")); if ((p_e [e.getIndex ()] < control_minLinkPersistenceProb.getDouble() - 1E-3) ||(p_e [e.getIndex ()] > 1E-3 + control_maxNodePersistenceProb.getDouble())) throw new RuntimeException ("Bad"); }
379                double [] q_n = new double [N]; 
380                for (Node n : currentNetPlan.getNodes ()) 
381                {
382                        for (Link e : n.getOutgoingLinks()) q_n [n.getIndex ()] += p_e [e.getIndex ()];
383                        if (Math.abs (q_n [n.getIndex ()] - Double.parseDouble (n.getAttribute("q_n"))) > 1E-3) throw new RuntimeException ("Bad");
384                        if (q_n [n.getIndex ()] > control_maxNodePersistenceProb.getDouble() + 1E-3) throw new RuntimeException ("Bad");
385                }
386                for (Link e : currentNetPlan.getLinks ())
387                {
388                        final double u_e = e.getCapacity();
389                        double supposedCapacity = control_linkNominalCapacities.get(e.getIndex ()) * control_p_e.get (e.getIndex ());
390                        for (Node n : control_nodesInterfTo_e.get(e)) supposedCapacity *= 1 - q_n [n.getIndex ()];
391                        if (Math.abs (u_e - supposedCapacity) > 1e-3) throw new RuntimeException ("Bad");
392                }
393        }
394        
395}