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.Arrays;
017import java.util.HashMap;
018import java.util.List;
019import java.util.Map;
020import java.util.Random;
021
022import cern.colt.matrix.tdouble.DoubleFactory1D;
023import cern.colt.matrix.tdouble.DoubleFactory2D;
024import cern.colt.matrix.tdouble.DoubleMatrix1D;
025import cern.colt.matrix.tdouble.DoubleMatrix2D;
026
027import com.jom.OptimizationProblem;
028import com.net2plan.interfaces.networkDesign.Demand;
029import com.net2plan.interfaces.networkDesign.Link;
030import com.net2plan.interfaces.networkDesign.Net2PlanException;
031import com.net2plan.interfaces.networkDesign.NetPlan;
032import com.net2plan.interfaces.networkDesign.Node;
033import com.net2plan.interfaces.networkDesign.Route;
034import com.net2plan.interfaces.simulation.IEventProcessor;
035import com.net2plan.interfaces.simulation.SimEvent;
036import com.net2plan.utils.Constants.RoutingType;
037import com.net2plan.utils.DoubleUtils;
038import com.net2plan.utils.GradientProjectionUtils;
039import com.net2plan.utils.InputParameter;
040import com.net2plan.utils.Quadruple;
041import com.net2plan.utils.TimeTrace;
042import com.net2plan.utils.Triple;
043
044/** 
045 * This module implements a distributed dual-gradient based algorithm, for iteratively adapting the network routing.
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_sec10_2_adaptiveRoutingDual.m">{@code fig_sec10_2_adaptiveRoutingDual.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 Flow assignment (FA), Distributed algorithm, Dual gradient algorithm
056 * @net2plan.ocnbooksections Section 10.2
057 * @net2plan.inputParameters 
058 * @author Pablo Pavon-Marino
059 */
060public class Online_evProc_adaptiveRoutingDual 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 InputParameter signaling_isSynchronous = new InputParameter ("signaling_isSynchronous", false , "true if all the distributed agents involved wake up synchronously to send the signaling messages");
066        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);
067        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);
068        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);
069        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);
070        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);
071        private InputParameter update_isSynchronous = new InputParameter ("update_isSynchronous", false , "true if all the distributed agents involved wake up synchronousely to update its state");
072        private InputParameter update_averageInterUpdateTime = new InputParameter ("update_averageInterUpdateTime", 1.0 , "Average time between two updates of an agent" , 0 , false , Double.MAX_VALUE , true);
073        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);
074        private InputParameter simulation_maxNumberOfUpdateIntervals = new InputParameter ("simulation_maxNumberOfUpdateIntervals", 500.0 , "Maximum number of update intervals in average per agent" , 0 , false , Double.MAX_VALUE , true);
075        private InputParameter simulation_randomSeed = new InputParameter ("simulation_randomSeed", (long) 1 , "Seed of the random number generator");
076        private InputParameter simulation_outFileNameRoot = new InputParameter ("simulation_outFileNameRoot", "adaptiveRoutingDual" , "Root of the file name to be used in the output files. If blank, no output");
077        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);
078        private InputParameter gradient_gammaStep = new InputParameter ("gradient_gammaStep", 0.001 , "Gamma step in the gradient algorithm" , 0 , false , Double.MAX_VALUE , true);
079        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);
080        private InputParameter gradient_maxGradientCoordinateChange = new InputParameter ("gradient_maxGradientCoordinateChange", 1.0e10 , "Maximum change in an iteration of a gradient coordinate" , 0 , false , Double.MAX_VALUE , true);
081        private InputParameter gradient_regularizationEpsilon = new InputParameter ("gradient_regularizationEpsilon", 1.0e-3 , "Regularization factor in the objective function term to make it strictly convex" , 0 , true , Double.MAX_VALUE , true);
082
083        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);
084        private InputParameter control_initialLinkWeights = new InputParameter ("control_initialLinkWeights", 0.5 , "Initial value of the link weights" , 0 , true , Double.MAX_VALUE , true);
085        
086        private Random rng;
087        private int E , R , N , D;
088
089        private static final int SIGNALING_WAKEUPTOSENDMESSAGE = 300;
090        private static final int SIGNALING_RECEIVEDMESSAGE = 301;
091        private static final int UPDATE_WAKEUPTOUPDATE = 302;
092
093        private NetPlan currentNetPlan;
094        private DoubleMatrix1D routing_price_e;
095        
096        private DoubleMatrix2D routing_mostUpdatedLinkPriceKnownByNode_ne; 
097        private DoubleMatrix1D routing_previousLinkWeight_e;
098        private int [][] control_routeIndexes_d;
099        
100        private TimeTrace stat_traceOf_xp;
101        private TimeTrace stat_traceOf_pie;
102        private TimeTrace stat_traceOf_ye;
103        private TimeTrace stat_traceOf_objFunction;
104
105        @Override
106        public String getDescription()
107        {
108                return "This module implements a distributed dual-gradient based algorithm, for iteratively adapting the network routing.";
109        }
110
111        @Override
112        public List<Triple<String, String, String>> getParameters()
113        {
114                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
115                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
116        }
117
118        @Override
119        public void initialize(NetPlan currentNetPlan, Map<String, String> algorithmParameters, Map<String, String> simulationParameters, Map<String, String> net2planParameters)
120        {
121                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
122                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
123
124                if (currentNetPlan.getNumberOfLayers() != 1) throw new Net2PlanException ("This algorithm works in single layer networks");
125
126                this.rng = new Random (simulation_randomSeed.getLong());
127                this.currentNetPlan = currentNetPlan;
128                this.E = currentNetPlan.getNumberOfLinks();
129                this.N = currentNetPlan.getNumberOfNodes();
130                this.D = currentNetPlan.getNumberOfDemands();
131                
132                /* Check all nodes have output links */
133                for (Node n : currentNetPlan.getNodes()) if (n.getOutgoingLinks().isEmpty()) throw new Net2PlanException ("All nodes should have output links");
134
135                /* Set the initial prices in the links to all zeros */
136                this.routing_price_e = DoubleFactory1D.dense.make (E , control_initialLinkWeights.getDouble());
137                this.routing_previousLinkWeight_e = routing_price_e.copy ();
138                for (Link e : currentNetPlan.getLinks ()) e.setAttribute("pi_e" , "" + routing_price_e.get (e.getIndex ()));
139                
140                /* Initialize the information each demand knows of the prices of all the links */
141                this.routing_mostUpdatedLinkPriceKnownByNode_ne = DoubleFactory2D.dense.make (N,E,control_initialLinkWeights.getDouble());
142                
143                /* Sets the initial routing, according to the prices, with all the traffic balanced equally in all the paths of the demand */
144                currentNetPlan.removeAllUnicastRoutingInformation();
145                currentNetPlan.setRoutingType(RoutingType.SOURCE_ROUTING);
146                this.currentNetPlan.addRoutesFromCandidatePathList(null , "K" , "" + control_maxNumberOfPathsPerDemand.getInt());
147                this.R = currentNetPlan.getNumberOfRoutes ();
148                this.control_routeIndexes_d = new int [D][];
149                for (Demand d : this.currentNetPlan.getDemands())
150                {
151                        control_routeIndexes_d [d.getIndex ()] = new int [d.getRoutes ().size ()]; 
152                        int counter_r = 0; for (Route r : d.getRoutes ()) control_routeIndexes_d [d.getIndex ()][counter_r ++] = r.getIndex ();
153                }
154
155                
156                for (Demand d : this.currentNetPlan.getDemands())
157                        updateDemandRoutingFromWeightsKnown (d);                
158
159                //System.out.println("Initial maximum link utilization: " + this.currentNetPlan.getLinkMaximumUtilization());
160                
161                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
162                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
163                /* Initially all nodes receive a "wake up to transmit" event, aligned at time zero or y asynchr => randomly chosen */
164                for (Link e : currentNetPlan.getLinks())
165                {
166                        final double signalingTime = (signaling_isSynchronous.getBoolean())? signaling_averageInterMessageTime.getDouble() : Math.max(0 , signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
167                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , e));
168                }
169                for (Node n : currentNetPlan.getNodes())
170                {
171                        final double updateTime = (update_isSynchronous.getBoolean())? update_averageInterUpdateTime.getDouble() : Math.max(0 , update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
172                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE , n));
173                }
174
175//              for (long e : currentNetPlan.getLinkIds())
176//              {
177//                      final double signalingTime = (routing_isSignalingSynchronous)? routing_averageSignalingInterval : Math.max(0 , routing_averageSignalingInterval + routing_maxFluctuationSignalingInterval * (rng.nextDouble() - 0.5));
178//                      scheduleEvent(new SimEvent (signalingTime , e , SimEvent.DestinationModule.EVENT_PROCESSOR).setEventType(ROUTING_LINK_UPDATEGRADIENTANDSIGNALTONODES));
179//              }
180//              for (long n : currentNetPlan.getNodeIds())
181//              {
182//                      final double recomputingTime = (routing_isRecomputingSynchronous)? routing_averageRecomputingInterval : Math.max(0 , routing_averageRecomputingInterval + routing_maxFluctuationRecomputingInterval * (rng.nextDouble() - 0.5));
183//                      scheduleEvent(new SimEvent (recomputingTime , n , SimEvent.DestinationModule.EVENT_PROCESSOR).setEventType(ROUTING_NODEWAKEUPTO_RECOMPUTEROUTING));
184//              }
185
186                /* Intialize the traces */
187                this.stat_traceOf_xp = new TimeTrace ();
188                this.stat_traceOf_pie = new TimeTrace ();
189                this.stat_traceOf_ye = new TimeTrace ();
190                this.stat_traceOf_objFunction = new TimeTrace (); 
191                this.stat_traceOf_xp.add(0.0, this.currentNetPlan.getVectorRouteCarriedTraffic());
192                this.stat_traceOf_pie.add(0.0, this.routing_price_e.copy ());
193                this.stat_traceOf_ye.add(0.0, this.currentNetPlan.getVectorLinkTotalCarriedTraffic());
194                this.stat_traceOf_objFunction.add(0.0, computeObjectiveFucntionFromNetPlan());
195
196                /* */
197//              System.out.println("U_e: " + this.currentNetPlan.getLinkCapacityMap());
198//              log.println("Initial rouetrs in this.currentNetPlan.getRouteCarriedTrafficMap(): " + this.currentNetPlan.getRouteCarriedTrafficMap());
199        }
200
201        @Override
202        public void processEvent(NetPlan currentNetPlan, SimEvent event)
203        {
204                final double t = event.getEventTime();
205                switch (event.getEventType())
206                {
207                case SIGNALING_RECEIVEDMESSAGE: // A node receives a price information from a link
208                {
209                        final Triple<Node,Link,Double> signalInfo = (Triple<Node,Link,Double>) event.getEventObject();
210                        final Node nMe = signalInfo.getFirst();
211                        final Link e = signalInfo.getSecond();
212                        final double link_price_e = signalInfo.getThird();
213                        this.routing_mostUpdatedLinkPriceKnownByNode_ne.set(nMe.getIndex () , e.getIndex () , link_price_e);
214                        break;
215                }
216                
217                case SIGNALING_WAKEUPTOSENDMESSAGE: // A link updates its price (gradient), and signals it to all the nodes
218                {
219                        final Link eMe = (Link) event.getEventObject();
220
221                        /* Update the gradient iteration */
222                        final double previous_pi_e = this.routing_previousLinkWeight_e.get(eMe.getIndex ());
223                        final double old_pi_e = this.routing_price_e.get(eMe.getIndex ());
224                        final double gradient_e = eMe.getCarriedTrafficIncludingProtectionSegments() - eMe.getCapacity() +  2*gradient_maxGradientAbsoluteNoise.getDouble()*(this.rng.nextDouble()-0.5);
225                        final double new_pi_e_notProjected = old_pi_e + this.gradient_gammaStep.getDouble() * gradient_e + this.gradient_heavyBallBetaParameter.getDouble() * (old_pi_e - previous_pi_e);
226                        double new_pi_e_projected = Math.max(0, new_pi_e_notProjected);
227                        if (gradient_maxGradientCoordinateChange.getDouble() > 0)
228                                new_pi_e_projected = GradientProjectionUtils.scaleDown_maxAbsoluteCoordinateChange (old_pi_e , new_pi_e_projected ,  gradient_maxGradientCoordinateChange.getDouble());
229
230//                      if (eIdMe == 0)
231//                      {
232//                              System.out.println("previous_pi_e: " + previous_pi_e);
233//                              System.out.println("old_pi_e: " + old_pi_e);
234//                              System.out.println("gradient_e: " + gradient_e);
235//                              System.out.println("new_pi_e_notProjected: " + new_pi_e_notProjected);
236//                              System.out.println("new_pi_e_projected: " + new_pi_e_projected);
237//                      }
238                                
239                        this.routing_previousLinkWeight_e.set(eMe.getIndex (), old_pi_e);
240                        this.routing_price_e.set(eMe.getIndex () , new_pi_e_projected);
241                        eMe.setAttribute("pi_e" , "" + routing_price_e.get (eMe.getIndex ()));
242
243                        
244                        /* Send the events of the signaling information messages to all the nodes */
245                        if (rng.nextDouble() >= this.signaling_signalingLossProbability.getDouble()) // the signaling may be lost => lost to all nodes
246                                for (Node n : currentNetPlan.getNodes ())
247                                {
248                                        final double signalingReceptionTime = t + Math.max(0 , signaling_averageDelay.getDouble() + signaling_maxFluctuationInDelay.getDouble() * (rng.nextDouble() - 0.5));
249                                        this.scheduleEvent(new SimEvent (signalingReceptionTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_RECEIVEDMESSAGE , Triple.of(n , eMe , routing_price_e.get(eMe.getIndex ()))));
250                                }
251                        
252                        /* Re-schedule when to wake up again */
253                        final double signalingTime = signaling_isSynchronous.getBoolean()? t + signaling_averageInterMessageTime.getDouble() : Math.max(t , t + signaling_averageInterMessageTime.getDouble() + signaling_maxFluctuationInterMessageTime.getDouble() * (rng.nextDouble() - 0.5));
254                        this.scheduleEvent(new SimEvent (signalingTime , SimEvent.DestinationModule.EVENT_PROCESSOR , SIGNALING_WAKEUPTOSENDMESSAGE , eMe));
255                        break;
256                        
257                }
258                
259                case UPDATE_WAKEUPTOUPDATE: // a node updates its source routing using the most updated information it has
260                {
261                        final Node nMe = (Node) event.getEventObject(); 
262                        for (Demand d : nMe.getOutgoingDemands())
263                                updateDemandRoutingFromWeightsKnown (d);                
264
265                        final double updateTime = update_isSynchronous.getBoolean()? t + update_averageInterUpdateTime.getDouble() : Math.max(t , t + update_averageInterUpdateTime.getDouble() + update_maxFluctuationInterUpdateTime.getDouble() * (rng.nextDouble() - 0.5));
266                        this.scheduleEvent(new SimEvent (updateTime , SimEvent.DestinationModule.EVENT_PROCESSOR , UPDATE_WAKEUPTOUPDATE,  nMe));
267
268                        this.stat_traceOf_xp.add(t, this.currentNetPlan.getVectorRouteCarriedTraffic());
269                        this.stat_traceOf_pie.add(t, this.routing_price_e.copy ());
270                        this.stat_traceOf_ye.add(t, this.currentNetPlan.getVectorLinkTotalCarriedTraffic());
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_pie.printToFile(new File (simulation_outFileNameRoot.getString() + "_pie.txt"));
290                stat_traceOf_ye.printToFile(new File (simulation_outFileNameRoot.getString() + "_ye.txt"));
291                stat_traceOf_objFunction.printToFile(new File (simulation_outFileNameRoot.getString() + "_objFunc.txt"));
292                Quadruple<DoubleMatrix1D,DoubleMatrix1D,DoubleMatrix1D,Double> pair = computeOptimumSolution (true);
293                DoubleMatrix1D x_p_opt = pair.getFirst();
294                DoubleMatrix1D pi_e = pair.getSecond();  
295                DoubleMatrix1D y_e = pair.getThird();  
296                double optCost = pair.getFourth();  
297                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jomReg_objFunc.txt"), optCost);
298                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jomReg_xp.txt"), x_p_opt);
299                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jomReg_pie.txt"), pi_e);
300                TimeTrace.printToFile(new File (simulation_outFileNameRoot.getString() + "_jomReg_ye.txt"), y_e);
301                return null;
302        }
303        
304        private double computeObjectiveFucntionFromNetPlan ()
305        {
306                double objFunc = 0;
307                for (Route r : this.currentNetPlan.getRoutes())
308                {
309                        final double x_r = r.getCarriedTraffic();
310                        objFunc +=  x_r * r.getNumberOfHops() + this.gradient_regularizationEpsilon.getDouble() * Math.pow(x_r, 2);
311                }
312                return objFunc;
313        }
314        
315        private Quadruple<DoubleMatrix1D,DoubleMatrix1D,DoubleMatrix1D,Double> computeOptimumSolution (boolean isRegularized)
316        {
317                /* Modify the map so that it is the pojection where all elements sum h_d, and are non-negative */
318                final int P = this.currentNetPlan.getNumberOfRoutes();
319                final int E = this.currentNetPlan.getNumberOfLinks();
320                final DoubleMatrix1D h_d = this.currentNetPlan.getVectorDemandOfferedTraffic();
321                final DoubleMatrix1D u_e = this.currentNetPlan.getVectorLinkCapacity();
322                final DoubleMatrix1D h_p = this.currentNetPlan.getVectorRouteOfferedTrafficOfAssociatedDemand();
323                                
324                OptimizationProblem op = new OptimizationProblem();
325                op.addDecisionVariable("x_p", false , new int[] { 1 , P }, DoubleUtils.zeros(P) , h_p.toArray()); 
326                op.addDecisionVariable("y_e", false , new int[] { 1 , E }, DoubleUtils.zeros(E) , u_e.toArray()); 
327                op.setInputParameter("h_d", h_d , "row");
328                op.setInputParameter("u_e", u_e , "row");
329                op.setInputParameter("epsilon", this.gradient_regularizationEpsilon.getDouble());
330                
331                /* Set some input parameters */
332                if (!isRegularized)
333                        op.setObjectiveFunction("minimize", "sum (y_e)");
334                else
335                        op.setObjectiveFunction("minimize", "sum (y_e) + epsilon * x_p * x_p'");
336
337                /* VECTORIAL FORM OF THE CONSTRAINTS  */
338                op.setInputParameter("A_dp", currentNetPlan.getMatrixDemand2RouteAssignment());
339                op.setInputParameter("A_ep", currentNetPlan.getMatrixLink2RouteAssignment());
340                op.addConstraint("A_dp * x_p' == h_d'"); // for each demand, the 100% of the traffic is carried (summing the associated paths)
341                op.addConstraint("A_ep * x_p' == y_e'"); // set the traffic in each link equal to y_e
342                
343                /* Call the solver to solve the problem */
344                op.solve("ipopt");
345
346                /* If an optimal solution was not found, quit */
347                if (!op.solutionIsOptimal ()) throw new RuntimeException ("An optimal solution was not found");
348        
349                /* Retrieve the optimum solutions. Convert the bps into Erlangs */
350                final DoubleMatrix1D x_p = op.getPrimalSolution("x_p").view1D();
351                final DoubleMatrix1D y_e = op.getPrimalSolution("y_e").view1D();
352                final DoubleMatrix1D pi_e = op.getMultiplierOfUpperBoundConstraintToPrimalVariables("y_e").view1D();
353                
354                /* Check the solution */
355                for (Link e : this.currentNetPlan.getLinks())
356                {
357                        double yThisLink = 0; for (Route r : e.getTraversingRoutes()) yThisLink += x_p.get(r.getIndex ());
358                        if (yThisLink > e.getCapacity() + 1E-3) throw new RuntimeException ("Bad");
359                }
360                for (Demand d : this.currentNetPlan.getDemands())
361                {
362                        double trafDemand = 0; for (Route r : d.getRoutes()) trafDemand += x_p.get(r.getIndex ());
363                        if (Math.abs(trafDemand - d.getOfferedTraffic()) > 1E-3) throw new RuntimeException ("Bad");
364                }
365                
366                /* Compute optimum cost */
367                double objFunc = 0; for (Route r : this.currentNetPlan.getRoutes()) 
368                {
369                        final double x_r = x_p.get(r.getIndex ()); 
370                        objFunc +=  x_r * r.getNumberOfHops();
371                        if (isRegularized) objFunc += this.gradient_regularizationEpsilon.getDouble() * Math.pow(x_r, 2);
372                }
373                return Quadruple.of(x_p,pi_e,y_e,objFunc);
374
375        }
376        
377        
378        private void updateDemandRoutingFromWeightsKnown (Demand d)
379        {
380                final Node a_d = d.getIngressNode();
381                final double h_d = d.getOfferedTraffic();
382                final DoubleMatrix1D weightsKnown_e = this.routing_mostUpdatedLinkPriceKnownByNode_ne.viewRow(a_d.getIndex ()); 
383                DoubleMatrix1D a_k = DoubleFactory1D.dense.make (R);
384//              Map<Long,Double> map_a_k = new HashMap<Long,Double> ();
385                for (Route r : d.getRoutes())
386                {
387                        double val = 0; 
388                        for (Link e : r.getSeqLinksRealPath()) val += weightsKnown_e.get(e.getIndex()) + 1; 
389                        a_k.set(r.getIndex (), val);
390//                      map_a_k.put (r.getId () , val);
391                }
392
393//              Map<Long,Double> x_r_map = GradientProjectionUtils.regularizedProjection_sumEquality (map_a_k , null , null , h_d , this.gradient_regularizationEpsilon.getDouble());
394//              for (Route r : d.getRoutes ())
395//                      r.setCarriedTraffic(x_r_map.get(r.getId ()) ,  x_r_map.get(r.getId ()));
396//              System.out.println ("MAP: x_r this route: " +  currentNetPlan.getVectorRouteCarriedTraffic().viewSelection(control_routeIndexes_d [d.getIndex ()])); 
397//
398//              System.out.println ("control_routeIndexes_d [d.getIndex ()]: " + Arrays.toString (control_routeIndexes_d [d.getIndex ()]));
399//              System.out.println ("a_k.viewSelection control_routeIndexes_d [d.getIndex ()]): " + a_k.viewSelection (control_routeIndexes_d [d.getIndex ()]));
400//              
401                DoubleMatrix1D x_r = GradientProjectionUtils.regularizedProjection_sumEquality (a_k , control_routeIndexes_d [d.getIndex ()] , null , h_d , this.gradient_regularizationEpsilon.getDouble());
402                for (Route r : d.getRoutes ())
403                        r.setCarriedTraffic(x_r.get(r.getIndex ()) ,  x_r.get(r.getIndex ()));
404
405//              System.out.println ("DOUBLE: x_r this route: " +  currentNetPlan.getVectorRouteCarriedTraffic().viewSelection(control_routeIndexes_d [d.getIndex ()])); 
406
407                
408//              if (1 == 1) throw new RuntimeException ("here");
409                
410//              if (d == this.currentNetPlan.getDemandIdsVector() [1])
411//              {
412//                      System.out.println("Demand " + d + ", h_d: " + h_d + ", epsilon: " + routing_regularizationEpsilon + ", a_k: : " + a_k);
413//                      System.out.println("Demand " + d + ", x_r: " + x_r);
414//                      System.out.println("y_e: " + this.currentNetPlan.getLinkCarriedTrafficMap());
415//                      System.out.println("pi_e: " + this.routing_price_e);
416//                      }
417
418        }
419        
420}