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}