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;
027
028import com.jom.OptimizationProblem;
029import com.net2plan.interfaces.networkDesign.Demand;
030import com.net2plan.interfaces.networkDesign.Link;
031import com.net2plan.interfaces.networkDesign.Net2PlanException;
032import com.net2plan.interfaces.networkDesign.NetPlan;
033import com.net2plan.interfaces.networkDesign.Node;
034import com.net2plan.interfaces.networkDesign.Route;
035import com.net2plan.interfaces.simulation.IEventProcessor;
036import com.net2plan.interfaces.simulation.SimEvent;
037import com.net2plan.utils.DoubleUtils;
038import com.net2plan.utils.GradientProjectionUtils;
039import com.net2plan.utils.InputParameter;
040import com.net2plan.utils.Pair;
041import com.net2plan.utils.TimeTrace;
042import com.net2plan.utils.Triple;
043import com.net2plan.utils.Constants.RoutingType;
044
045/** 
046 * This module implements a distributed primal-gradient based algorithm, for iteratively adapting the network routing.
047 *
048 * Ths event processor is adapted to permit observing the algorithm performances under user-defined conditions, 
049 * including asynchronous distributed executions, where signaling can be affected by losses and/or delays, and/or measurement errors. 
050 * The time evolution of different metrics can be stored in output files, for later processing. 
051 * As an example, see the <a href="../../../../../../graphGeneratorFiles/fig_sec9_3_adaptiveRoutingPrimal.m">{@code fig_sec9_3_adaptiveRoutingPrimal.m}</a> MATLAB file used for generating the graph/s of the case study in the 
052 * <a href="http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119013356.html">book</a> using this algorithm.
053 * 
054 * To simulate a network with this module, use the {@code Online_evGen_doNothing} generator.
055 * 
056 * @net2plan.keywords Flow assignment (FA), Distributed algorithm, Primal gradient algorithm
057 * @net2plan.ocnbooksections Section 9.3
058 * @net2plan.inputParameters 
059 * @author Pablo Pavon-Marino
060 */
061public class Online_evProc_adaptiveRoutingPrimal extends IEventProcessor 
062{
063        private InputParameter signaling_isSynchronous = new InputParameter ("signaling_isSynchronous", false , "true if all the distributed agents involved wake up synchronously to send the signaling messages");
064        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);
065        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);
066        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);
067        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);
068        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);
069        private InputParameter update_isSynchronous = new InputParameter ("update_isSynchronous", false , "true if all the distributed agents involved wake up synchronousely to update its state");
070        private InputParameter update_averageInterUpdateTime = new InputParameter ("update_averageInterUpdateTime", 1.0 , "Average time between two updates of an agent" , 0 , false , Double.MAX_VALUE , true);
071        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);
072        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);
073        private InputParameter gradient_gammaStep = new InputParameter ("gradient_gammaStep", 1.0 , "Gamma step in the gradient algorithm" , 0 , false , Double.MAX_VALUE , true);
074        private InputParameter gradient_heavyBallBetaParameter = new InputParameter ("gradient_heavyBallBetaParameter", 0.0 , "Beta parameter of heavy ball, between 0 and 1. Value 0 means no heavy ball" , 0 , true , 1.0 , true);
075        private InputParameter gradient_maxGradientCoordinateChange = new InputParameter ("gradient_maxGradientCoordinateChange", 1000.0 , "Maximum change in an iteration of a gradient coordinate" , 0 , false , Double.MAX_VALUE , true);
076        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);
077        private InputParameter simulation_randomSeed = new InputParameter ("simulation_randomSeed", (long) 1 , "Seed of the random number generator");
078        private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "adaptiveRoutingPrimal" , "Root of the file name to be used in the output files. If blank, no output");
079        private InputParameter control_maxNumberOfPathsPerDemand = new InputParameter ("control_maxNumberOfPathsPerDemand", (int) 5 , "Number of acceptable paths per demand. Take the k-shortest paths in number of hops" , 1 , Integer.MAX_VALUE);
080
081        private static PrintStream getNulFile () { try { return new PrintStream (new FileOutputStream ("NUL") , false); } catch (Exception e) {e.printStackTrace(); throw new RuntimeException ("Not NUL file"); }   } 
082        private PrintStream log = getNulFile (); //System.err;//new PrintStream (new FileOutputStream ("NUL") , true);
083
084        private Random rng;
085
086        private static final int SIGNALING_WAKEUPTOSENDMESSAGE = 300;
087        private static final int SIGNALING_RECEIVEDMESSAGE = 301;
088        private static final int UPDATE_WAKEUPTOUPDATE = 302;
089
090        private NetPlan currentNetPlan;
091
092        private DoubleMatrix1D control_price_e;
093        private DoubleMatrix2D control_mostUpdated_price_e_knownByNode_n;
094        private DoubleMatrix1D control_previousXp;
095        private int [][] control_routeIndexes_d;
096        private int E , R , N , D;
097        
098//      private Map<Node,DoubleMatrix1D> routing_mostUpdated_price_e_knownByNode_n; 
099//      private Map<Long,Set<Long>> routing_outgoingDemands_n;
100//      private Map<Long,Double> routing_previousXp;
101        
102        private TimeTrace stat_traceOf_xp;
103        private TimeTrace stat_traceOf_objFunction;
104        
105        
106        @Override
107        public String getDescription()
108        {
109                return "This module implements a distributed primal-gradient based algorithm, for iteratively adapting the network routing.";
110        }
111
112        @Override
113        public List<Triple<String, String, String>> getParameters()
114        {
115                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
116                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
117        }
118
119        @Override
120        public void initialize(NetPlan currentNetPlan, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
121        {
122                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
123                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
124
125                if (currentNetPlan.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks");
126
127                this.rng = new Random (simulation_randomSeed.getLong());
128                this.currentNetPlan = currentNetPlan;
129                this.E = currentNetPlan.getNumberOfLinks();
130                this.N = currentNetPlan.getNumberOfNodes();
131                this.D = currentNetPlan.getNumberOfDemands();
132                
133                /* Check all nodes have output links */
134                for (Node n : currentNetPlan.getNodes()) if (n.getOutgoingLinks().isEmpty()) throw new Net2PlanException ("All nodes should have output links");
135
136//              /* Compute the CPL */
137//              this.cpl = new CandidatePathList (this.currentNetPlan , "K" , "" + routing_maxNumberOfPathsPerDemand);
138                
139                /* Sets the initial routing, with all the traffic balanced equally in all the paths of the demand */
140                currentNetPlan.removeAllUnicastRoutingInformation();
141                currentNetPlan.setRoutingType(RoutingType.SOURCE_ROUTING);
142                this.currentNetPlan.addRoutesFromCandidatePathList(null , "K" , "" + control_maxNumberOfPathsPerDemand.getInt());
143                this.R = currentNetPlan.getNumberOfRoutes ();
144                this.control_previousXp = DoubleFactory1D.dense.make (currentNetPlan.getNumberOfRoutes());
145                this.control_routeIndexes_d = new int [D][];
146                for (Demand d : this.currentNetPlan.getDemands())
147                {
148                        final double h_d = d.getOfferedTraffic();
149                        final Set<Route> routes = d.getRoutes();
150                        if (routes.isEmpty()) throw new Net2PlanException ("A demand has no paths assigned: not connected topology");
151                        control_routeIndexes_d [d.getIndex ()] = new int [routes.size ()]; 
152                        int counter_r = 0;
153                        for (Route r : routes)
154                        {
155                                control_routeIndexes_d [d.getIndex ()][counter_r ++] = r.getIndex ();
156                                r.setCarriedTraffic(h_d / routes.size() , h_d / routes.size());
157                                control_previousXp.set(r.getIndex(), h_d / routes.size());
158                                log.println("Initial route: " + r + " demand: " +d + ", h_d: " + h_d + ", h_r = " + r.getCarriedTraffic() + ", seqLinks: "+ r.getSeqLinksRealPath());
159                        }
160                }
161
162                /* Set the initial prices in the nodes */
163                this.control_price_e = DoubleFactory1D.dense.make (E);
164                for (Link e : currentNetPlan.getLinks ()) this.control_price_e.set (e.getIndex() , computeLinkPriceFromNetPlan (e));
165                
166                /* Initialize the information each demand knows of the prices of all the links */
167                this.control_mostUpdated_price_e_knownByNode_n = DoubleFactory2D.dense.make (N , E);
168                for (int n = 0; n < N ; n ++) control_mostUpdated_price_e_knownByNode_n.viewRow (n).assign (control_price_e);
169                
170                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
171                for (Node n : currentNetPlan.getNodes())
172                {
173                        final double signalingTime = (signaling_isSynchronous.getBoolean())? signaling_averageInterMessageTime.getDouble() : Math.max(0 , signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
174                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , n));
175                        final double updateTime = (update_isSynchronous.getBoolean())? update_averageInterUpdateTime.getDouble() : Math.max(0 , update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
176                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE , n));
177                }
178
179                /* Intialize the traces */
180                this.stat_traceOf_xp = new TimeTrace();
181                this.stat_traceOf_objFunction = new TimeTrace (); 
182                this.stat_traceOf_xp.add(0.0, this.currentNetPlan.getVectorRouteCarriedTraffic());
183                this.stat_traceOf_objFunction.add(0.0, computeObjectiveFucntionFromNetPlan());
184
185                /* */
186                log.println("Initial rouetrs in this.currentNetPlan.getRouteCarriedTrafficMap(): " + this.currentNetPlan.getVectorRouteCarriedTraffic());
187        }
188
189        @Override
190        public void processEvent(NetPlan currentNetPlan, SimEvent event)
191        {
192                final double t = event.getEventTime();
193                switch (event.getEventType())
194                {
195                case SIGNALING_RECEIVEDMESSAGE: // A node receives from an out neighbor the q_nt for any destination
196                {
197                        final Pair<Node,Map<Link,Double>> signalInfo = (Pair<Node,Map<Link,Double>>) event.getEventObject();
198                        final Node nMe = signalInfo.getFirst();
199                        final Map<Link,Double> receivedInfo_price_e = signalInfo.getSecond();
200                        for (Link e : receivedInfo_price_e.keySet())
201                                this.control_mostUpdated_price_e_knownByNode_n.set (nMe.getIndex () , e.getIndex() , receivedInfo_price_e.get(e));
202                        log.println("t=" + t + ": SIGNALING_RECEIVEDMESSAGE node " + nMe);
203                        break;
204                }
205                
206                case SIGNALING_WAKEUPTOSENDMESSAGE: // A node broadcasts signaling info to its 1 hop neighbors
207                {
208                        final Node nMe = (Node) event.getEventObject();
209                        log.println("t=" + t + ": ROUTING_NODEWAKEUPTO_TRANSMITSIGNAL node " + nMe);
210
211                        /* Create the info I will signal */
212                        Map<Link,Double> infoToSignal_price_e = new HashMap<Link,Double> ();
213                        for (Link e : nMe.getOutgoingLinks())
214                                infoToSignal_price_e.put (e , this.computeLinkPriceFromNetPlan(e));
215
216                        /* Send the events of the signaling information messages to all the nodes */
217                        if (rng.nextDouble() >= this.signaling_signalingLossProbability.getDouble()) // the signaling may be lost => lost to all nodes
218                                for (Node n : currentNetPlan.getNodes ())
219                                {
220                                        final double signalingReceptionTime = t + Math.max(0 , signaling_averageDelay.getDouble() + signaling_maxFluctuationInDelay.getDouble() * (rng.nextDouble() - 0.5));
221                                        this.scheduleEvent(new SimEvent (signalingReceptionTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_RECEIVEDMESSAGE , Pair.of(n , infoToSignal_price_e)));
222                                }
223                        
224                        /* Re-schedule when to wake up again */
225                        final double signalingTime = signaling_isSynchronous.getBoolean()? t + signaling_averageInterMessageTime.getDouble() : Math.max(t , t + signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
226                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , nMe));
227                        break;
228                }
229
230                case UPDATE_WAKEUPTOUPDATE: // a node updates its p_n, p_e values, using most updated info available
231                {
232                        final Node nMe = (Node) event.getEventObject(); 
233                        log.println("t=" + t + ": ROUTING_NODEWAKEUPTO_RECOMPUTE node " + nMe);
234                        //log.println("p_e values Before = " + mac_p_e);
235                        
236                        //if (1 ==1) break;
237                        DoubleMatrix1D x_p = this.currentNetPlan.getVectorRouteCarriedTraffic();
238                        DoubleMatrix1D infoIKnow_price_e = this.control_mostUpdated_price_e_knownByNode_n.viewRow (nMe.getIndex ());
239                        for (Demand d : nMe.getOutgoingDemands())
240                        {
241                                final double h_d = d.getOfferedTraffic();
242                                DoubleMatrix1D new_x_p = DoubleFactory1D.dense.make (R);
243                                for (int rIndex : control_routeIndexes_d [d.getIndex ()])
244                                {
245                                        final Route r = currentNetPlan.getRoute (rIndex);
246                                        final double old_xp = x_p.get(rIndex);
247                                        double accumPrices = 0;
248                                        for (Link e : r.getSeqLinksRealPath()) accumPrices += infoIKnow_price_e.get(e.getIndex());
249                                        accumPrices += 2*gradient_maxGradientAbsoluteNoise.getDouble()*(this.rng.nextDouble()-0.5);
250                                        new_x_p.set (r.getIndex () , old_xp - this.gradient_gammaStep.getDouble() * accumPrices + this.gradient_heavyBallBetaParameter.getDouble() * (r.getCarriedTraffic() - this.control_previousXp.get(r.getIndex ())  ));
251                                }
252
253                                GradientProjectionUtils.euclideanProjection_sumEquality (new_x_p , control_routeIndexes_d [d.getIndex ()] , h_d);
254
255                                if (gradient_maxGradientCoordinateChange.getDouble() > 0)
256                                        GradientProjectionUtils.scaleDown_maxAbsoluteCoordinateChange (x_p , new_x_p , control_routeIndexes_d [d.getIndex ()] ,  gradient_maxGradientCoordinateChange.getDouble());
257
258                                for (int rIndex : control_routeIndexes_d [d.getIndex ()])
259                                {
260                                        final double newCarriedTraffic = new_x_p.get (rIndex);
261                                        final Route r = currentNetPlan.getRoute (rIndex);
262                                        control_previousXp.set (rIndex , r.getCarriedTraffic());
263                                        r.setCarriedTraffic(newCarriedTraffic , newCarriedTraffic);
264                                }
265                        }
266
267                        final double updateTime = update_isSynchronous.getBoolean()? t + update_averageInterUpdateTime.getDouble() : Math.max(t , t + update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
268                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE,  nMe));
269
270                        this.stat_traceOf_xp.add(t, currentNetPlan.getVectorRouteCarriedTraffic());
271                        this.stat_traceOf_objFunction.add(t, computeObjectiveFucntionFromNetPlan());
272
273                        if (t > this.simulation_maxNumberOfUpdateIntervals.getDouble() * this.update_averageInterUpdateTime.getDouble()) { this.endSimulation (); }
274                        
275                        break;
276                }
277                        
278
279                default: throw new RuntimeException ("Unexpected received event");
280                }
281                
282                
283        }
284
285        public String finish (StringBuilder st , double simTime)
286        {
287                if (simulation_outFileNameRoot.getString().equals("")) return null;
288                stat_traceOf_xp.printToFile(new File (simulation_outFileNameRoot.getString() + "_xp.txt"));
289                stat_traceOf_objFunction.printToFile(new File (simulation_outFileNameRoot.getString() + "_objFunc.txt"));
290                Pair<DoubleMatrix1D,Double> optPair = computeOptimumSolution ();
291                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_xp.txt"), optPair.getFirst());
292                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jom_objFunc.txt"), optPair.getSecond());
293                return null;
294        }
295        
296        private double computeLinkPriceFromNetPlan (Link e)
297        {
298                final double y_e = e.getCarriedTrafficIncludingProtectionSegments();
299                final double u_e = e.getCapacity();
300                if (u_e == 0) throw new RuntimeException ("Zero capacity in a link means");
301                final double gradient = (y_e / u_e > 0.99)? 1/Math.pow(u_e * 0.01,2) : 1/Math.pow(u_e - y_e,2); 
302                return gradient;
303        }
304
305        
306        
307        private double computeObjectiveFucntionFromNetPlan ()
308        {
309                double objFunc = 0;
310                for (Link e : this.currentNetPlan.getLinks())
311                {
312                        final double y_e = e.getCarriedTrafficIncludingProtectionSegments();
313                        final double u_e = e.getCapacity();
314                        if (y_e / u_e > 0.99)
315                                objFunc += 1/(0.01*u_e) + 1/Math.pow(u_e * 0.01,2) * (y_e - 0.99*u_e);
316                        else
317                                objFunc += 1/(u_e - y_e);
318                }
319                return objFunc;
320        }
321        
322        private Pair<DoubleMatrix1D,Double> computeOptimumSolution ()
323        {
324                /* Modify the map so that it is the pojection where all elements sum h_d, and are non-negative */
325                final int R = this.currentNetPlan.getNumberOfRoutes();
326                final int E = this.currentNetPlan.getNumberOfLinks();
327                final DoubleMatrix1D h_d = this.currentNetPlan.getVectorDemandOfferedTraffic();
328                final DoubleMatrix1D u_e = this.currentNetPlan.getVectorLinkCapacity();
329                final DoubleMatrix1D h_p = this.currentNetPlan.getVectorRouteOfferedTrafficOfAssociatedDemand();
330                                
331                OptimizationProblem op = new OptimizationProblem();
332                op.addDecisionVariable("x_p", false , new int[] { 1 , R }, DoubleUtils.zeros(R) , h_p.toArray());
333                op.addDecisionVariable("y_e", false , new int[] { 1 , E }, DoubleUtils.zeros(E) , u_e.toArray()); 
334                op.setInputParameter("h_d", h_d , "row");
335                op.setInputParameter("h_p", h_p , "row");
336                op.setInputParameter("u_e", u_e , "row");
337                
338                /* Set some input parameters */
339                op.setObjectiveFunction("minimize", "sum (1./(u_e - y_e))");
340
341                /* VECTORIAL FORM OF THE CONSTRAINTS  */
342                op.setInputParameter("A_dp", currentNetPlan.getMatrixDemand2RouteAssignment());
343                op.setInputParameter("A_ep", currentNetPlan.getMatrixLink2RouteAssignment());
344                op.addConstraint("A_dp * x_p' == h_d'"); // for each demand, the 100% of the traffic is carried (summing the associated paths)
345                op.addConstraint("A_ep * x_p' == y_e'"); // set the traffic in each link equal to y_e
346                
347                /* Call the solver to solve the problem */
348                op.solve("ipopt");
349
350                /* If an optimal solution was not found, quit */
351                if (!op.solutionIsOptimal ()) throw new RuntimeException ("An optimal solution was not found");
352        
353                /* Retrieve the optimum solutions. Convert the bps into Erlangs */
354                final double optCost = op.getOptimalCost();
355                DoubleMatrix1D x_p = op.getPrimalSolution("x_p").view1D();
356//              Map<Long,Double> x_p_map = new HashMap<Long,Double> ();
357//              for (long p : cpl.getPathIds()) x_p_map.put(p, x_p_array [id2IndexPathMap.get(p)]); 
358
359                /* Check the solution */
360                for (Link e : this.currentNetPlan.getLinks())
361                {
362                        double y_e = 0; for (Route r : e.getTraversingRoutes()) y_e += x_p.get (r.getIndex());
363                        if (y_e > e.getCapacity() + 1E-3) throw new RuntimeException ("Bad");
364                }
365                for (Demand d : this.currentNetPlan.getDemands())
366                {
367                        double trafDemand = 0; for (Route r : d.getRoutes()) trafDemand += x_p.get(r.getIndex ());
368                        if (Math.abs(trafDemand - d.getOfferedTraffic()) > 1E-3) throw new RuntimeException ("Bad");
369                }
370                
371                return Pair.of(x_p,optCost);
372        }
373        
374
375        
376}