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 ******************************************************************************/
008
009
010
011
012package com.net2plan.examples.ocnbook.offline;
013
014import java.io.File;
015import java.io.PrintWriter;
016import java.util.ArrayList;
017import java.util.Collection;
018import java.util.Collections;
019import java.util.HashMap;
020import java.util.HashSet;
021import java.util.LinkedList;
022import java.util.List;
023import java.util.Map;
024import java.util.Map.Entry;
025import java.util.Random;
026import java.util.Set;
027
028import cern.colt.matrix.tdouble.DoubleFactory1D;
029import cern.colt.matrix.tdouble.DoubleFactory2D;
030import cern.colt.matrix.tdouble.DoubleMatrix1D;
031import cern.colt.matrix.tdouble.DoubleMatrix2D;
032import cern.jet.math.tdouble.DoubleFunctions;
033
034import com.jom.DoubleMatrixND;
035import com.jom.OptimizationProblem;
036import com.net2plan.interfaces.networkDesign.Demand;
037import com.net2plan.interfaces.networkDesign.IAlgorithm;
038import com.net2plan.interfaces.networkDesign.Link;
039import com.net2plan.interfaces.networkDesign.Net2PlanException;
040import com.net2plan.interfaces.networkDesign.NetPlan;
041import com.net2plan.interfaces.networkDesign.NetworkElement;
042import com.net2plan.interfaces.networkDesign.Node;
043import com.net2plan.interfaces.networkDesign.ProtectionSegment;
044import com.net2plan.interfaces.networkDesign.Route;
045import com.net2plan.interfaces.networkDesign.SharedRiskGroup;
046import com.net2plan.utils.Constants.RoutingType;
047import com.net2plan.utils.Constants.SearchType;
048import com.net2plan.utils.DoubleUtils;
049import com.net2plan.utils.InputParameter;
050import com.net2plan.utils.IntUtils;
051import com.net2plan.utils.Pair;
052import com.net2plan.utils.Quadruple;
053import com.net2plan.utils.Triple;
054
055/**
056 * This algorithm is devoted to solve the several network planning problems in an optical WDM network (fiber placement, RWA, under different recovery schemes), appearing in the case study in the book section mentioned below. 
057 * @net2plan.description 
058 * @net2plan.keywords WDM, Topology assignment (TA), Flow assignment (FA), GRASP, JOM
059 * @net2plan.ocnbooksections Section 12.10
060 * @net2plan.inputParameters 
061 * @author Pablo Pavon-Marino
062 */
063public class Offline_tcfa_wdmPhysicalDesign_graspAndILP implements IAlgorithm
064{
065        /* Stat */
066        private double stat_totalCost , stat_costinks , stat_costCirc , stat_costNodes;
067        private int stat_numLinks , stat_numCirc , stat_numNodesDeg2 , stat_numNodesDeg3;
068        private int stat_numItGRASP , stat_numLSIterationsReducingCost , stat_numSolverCalls , stat_numSolverCallsOk;
069        private double stat_totalTimeSecs , stat_totalTimeCreateProgramsSecs , stat_totalTimeSolverSecs;
070        
071        /* Parameters */
072        private DoubleMatrix1D linkCost_e;
073        private Random rng;
074        private NetPlan netPlan;
075        
076        /* Control */
077        private int Efm, N, nSRGs , D , R; 
078        private DoubleMatrix1D numCircH_d; 
079        private DoubleMatrix2D A_er , A_dr , A_rs , A_se , Aout_ne , Ain_ne, Abid_dd , Abid_rr , Abid_ee;
080        private Map<Demand,Demand> opposite_d;
081        private Map<Route,Route> opposite_r;
082        private int [][] routeList_d; // for each demand, the array of route ids assigned to it
083        private int [][] traversedLinks_r; // for each demand, the array of route ids assigned to it
084
085        
086        private InputParameter algorithm_randomSeed = new InputParameter ("algorithm_randomSeed", (long) 1 , "Seed of the random number generator");
087        private InputParameter algorithm_outputFileNameRoot = new InputParameter ("algorithm_outputFileNameRoot", "tcfaWDM_graspAndILP" , "Root of the file name to be used in the output files. If blank, no output");
088        private InputParameter algorithm_maxExecutionTimeInSeconds = new InputParameter ("algorithm_maxExecutionTimeInSeconds", (double) 300 , "Algorithm maximum running time in seconds" , 0 , false , Double.MAX_VALUE , true);
089        private InputParameter tcfa_circuitCost = new InputParameter ("tcfa_circuitCost", (double) 3.0 , "Cost of each circuit" , 0 , true , Double.MAX_VALUE , true);
090        private InputParameter tcfa_circuitCapacity_Gbps = new InputParameter ("tcfa_circuitCapacity_Gbps", (double) 10.0 , "Capacity of an optical circuit" , 0 , true , Double.MAX_VALUE , true);
091        private InputParameter tcfa_srgType = new InputParameter ("tcfa_srgType" , "#select# perBidirectionalLinkBundle noFailure", "Determines how the SRGs are initialized. The design must be tolerant to single SRG failures");
092        private InputParameter tcfa_srgMttfPer1000Km_hrs = new InputParameter ("tcfa_srgMttfPer1000Km_hrs" , (double) 4380 , "Mean Time To Fail (MTTF) of the SRGs defined is this value multiplied by the link length divided by 1000",  0 , false , Double.MAX_VALUE , true);
093        private InputParameter tcfa_srgMttr_hrs = new InputParameter ("tcfa_srgMttr_hrs" , (double) 24 , "Mean Time To Repair (MTTR) of the SRGs",  0 , true , Double.MAX_VALUE , true);
094        private InputParameter tcfa_recoveryType = new InputParameter ("tcfa_recoveryType" , "#select# 1+1 shared restoration", "Determines the type of network recovery mechanism in the network.");
095        private InputParameter tcfa_maxPathLengthInKm = new InputParameter ("tcfa_maxPathLengthInKm" , (double) 10000 , "Maximum length that a lightpath can have",  0 , true , Double.MAX_VALUE , true);
096        private InputParameter tcfa_maxPathNumberOfHops = new InputParameter ("tcfa_maxPathNumberOfHops" , (int) 4 , "Maximum number of hops that a lightpath can make",  1 , Integer.MAX_VALUE);
097        private InputParameter tcfa_maxNumberPathsPerDemand = new InputParameter ("tcfa_maxNumberPathsPerDemand" , (int) 200 , "Maximum number of paths in the candidate path list, for each demand",  1 , Integer.MAX_VALUE);
098        private InputParameter tcfa_linkCapacity_numCirc = new InputParameter ("tcfa_linkCapacity_numCirc" , (int) 80 , "Number of wavelengths available per fiber",  1 , Integer.MAX_VALUE);
099        private InputParameter tcfa_linkCostPerKm = new InputParameter ("tcfa_linkCostPerKm" , (double) 0.4 , "Cost of one km of optical fiber",  0 , true  , Double.MAX_VALUE , true);
100        private InputParameter tcfa_costNodeDegree2 = new InputParameter ("tcfa_costNodeDegree2" , (double) 4.0 , "Cost of one OADM of degree lower or equal than two",  0 , true  , Double.MAX_VALUE , true);
101        private InputParameter tcfa_costNodeDegreeHigherThan2 = new InputParameter ("tcfa_costNodeDegreeHigherThan2" , (double) 4.0 , "Cost of one OADM of degree higher than two",  0 , true  , Double.MAX_VALUE , true);
102        private InputParameter tcfa_networkDiameterKmToNormalize = new InputParameter ("tcfa_networkDiameterKmToNormalize" , (double) 4500 , "The link lengths are normalized so the longest link has this value",  0 , true  , Double.MAX_VALUE , true);
103        private InputParameter tcfa_maxNumIterations = new InputParameter ("tcfa_maxNumIterations", (int) 10000 , "Maximum number of iterations" , 1 , Integer.MAX_VALUE);
104        private InputParameter algorithm_solverName = new InputParameter ("solverName", "#select# glpk cplex ipopt", "The solver name to be used by JOM. GLPK and IPOPT are free, CPLEX commercial. GLPK and CPLEX solve linear problems w/w.o integer contraints. IPOPT is can solve nonlinear problems (if convex, returns global optimum), but cannot handle integer constraints");
105        private InputParameter algorithm_solverLibraryName = new InputParameter ("solverLibraryName", "", "The solver library full or relative path, to be used by JOM. Leave blank to use JOM default.");
106        private InputParameter algorithm_maxSolverTimeInSeconds = new InputParameter ("maxSolverTimeInSeconds", (double) 2 , "Maximum time granted to the solver in to solve each problem instance. If this time expires, the solver returns the best solution found so far (if a feasible solution is found)");
107        private InputParameter tcfa_bidirectionalLinks = new InputParameter ("tcfa_bidirectionalLinks", true , "If true, the topology of fibers deployed is constrained to be bidirectional (the same number of fibers in both directions)");
108
109        @Override
110        public String executeAlgorithm(NetPlan netPlan, Map<String, String> algorithmParameters, Map<String, String> net2planParameters)
111        {
112                /* Initialize all InputParameter objects defined in this object (this uses Java reflection) */
113                InputParameter.initializeAllInputParameterFieldsOfObject(this, algorithmParameters);
114                if (netPlan.getNumberOfNodes() == 0) throw new Net2PlanException ("The input design has no nodes");
115                
116                
117                try{ 
118                this.netPlan = netPlan;
119                /* Initializations */
120
121                netPlan.removeAllLinks(); netPlan.setRoutingType(RoutingType.SOURCE_ROUTING);
122                this.N = netPlan.getNumberOfNodes();
123                this.rng = new Random (algorithm_randomSeed.getLong ());
124                
125                /* Initialize links */
126                for (Node n1 : netPlan.getNodes()) for (Node n2 : netPlan.getNodes()) if (n1 != n2) netPlan.addLink(n1, n2, 0.01, netPlan.getNodePairEuclideanDistance(n1, n2), 200000 , null);
127                final double linkLengthKmFromEuclideanDistanceFactor = tcfa_networkDiameterKmToNormalize.getDouble() / netPlan.getVectorLinkLengthInKm().getMaxLocation() [0];
128                for (Link e : netPlan.getLinks()) e.setLengthInKm(e.getLengthInKm() * linkLengthKmFromEuclideanDistanceFactor);
129                this.Efm = netPlan.getNumberOfLinks(); if (Efm != N*(N-1)) throw new RuntimeException ("Bad, Efm" + Efm + ", N: " + N);
130                this.linkCost_e = DoubleFactory1D.dense.make(Efm); for (Link e : netPlan.getLinks ()) linkCost_e.set(e.getIndex () , tcfa_linkCostPerKm.getDouble() * e.getLengthInKm());
131                
132                /* If there is no traffic, create uniform traffic matrix */
133                if (netPlan.getNumberOfDemands() == 0)
134                        for (Node n1 : netPlan.getNodes()) for (Node n2 : netPlan.getNodes()) if (n1 != n2) netPlan.addDemand(n1, n2, tcfa_circuitCapacity_Gbps.getDouble(), null);
135
136                
137                /* Convert the demand offered traffic in the upper multiple of circuit capacity */
138                DoubleMatrix2D trafficMatrix = netPlan.getMatrixNode2NodeOfferedTraffic();
139                netPlan.removeAllDemands();
140                for (Node n1 : netPlan.getNodes()) for (Node n2 : netPlan.getNodes()) if (n1.getIndex () > n2.getIndex ())
141                {
142                        final double maxTrafficBidir = Math.max(trafficMatrix.get(n1.getIndex (), n2.getIndex ()), trafficMatrix.get(n2.getIndex (), n1.getIndex ()));
143                        netPlan.addDemand(n1, n2, tcfa_circuitCapacity_Gbps.getDouble() * Math.ceil(maxTrafficBidir / tcfa_circuitCapacity_Gbps.getDouble()) , null); 
144                        netPlan.addDemand(n2, n1, tcfa_circuitCapacity_Gbps.getDouble() * Math.ceil(maxTrafficBidir / tcfa_circuitCapacity_Gbps.getDouble()) , null); 
145                }
146                        
147                /* If 1+1 then decouple demands in one per channel */
148                if (tcfa_recoveryType.getString ().equals("1+1"))
149                {
150                        NetPlan npcopy = netPlan.copy();
151                        netPlan.removeAllDemands();
152                        for (Demand d : npcopy.getDemands())
153                        {
154                                final int numChannels = (int) Math.round (d.getOfferedTraffic() / tcfa_circuitCapacity_Gbps.getDouble ());
155                                for (int cont = 0 ; cont < numChannels ; cont ++) netPlan.addDemand(netPlan.getNodeFromId (d.getIngressNode().getId ()), netPlan.getNodeFromId (d.getEgressNode().getId ()), tcfa_circuitCapacity_Gbps.getDouble () , d.getAttributes());
156                        }
157                }
158
159                /* Initialize routes and demands to make it bidirectional */
160                Pair<Map<Demand,Demand>,Map<Route,Route>> pair = initializeNetPlanLinksBidirDemandsAndRoutes (netPlan);
161                this.opposite_d = pair.getFirst();
162                this.opposite_r = pair.getSecond();
163
164                /* Initialize demand info */
165                this.D = netPlan.getNumberOfDemands();
166                this.numCircH_d = netPlan.getVectorDemandOfferedTraffic().assign(DoubleFunctions.div(this.tcfa_circuitCapacity_Gbps.getDouble ())).assign(DoubleFunctions.rint);
167                System.out.println("Total number of circuits to establish: " + numCircH_d.zSum());
168                
169                //if (1==1) throw new RuntimeException ("Stop");
170                this.R = netPlan.getNumberOfRoutes();
171                //System.out.println ("numCircH_d: " + numCircH_d);
172
173                /* Initialize SRG information */
174                initializeSRGs(netPlan);
175                this.nSRGs = netPlan.getNumberOfSRGs();
176                this.routeList_d = new int [D][]; for (Demand d : netPlan.getDemands ()) routeList_d [d.getIndex ()] = getIndexes (d.getRoutes ());
177                this.traversedLinks_r = new int [R][]; for (Route r : netPlan.getRoutes ()) traversedLinks_r [r.getIndex ()] = getIndexes (r.getSeqLinksRealPath());
178                
179                
180                /* Initialize aux arrays for speed-up */
181                this.A_dr = netPlan.getMatrixDemand2RouteAssignment();
182                for (Demand d : netPlan.getDemands ()) if (d.getRoutes().isEmpty()) throw new Net2PlanException ("A demand has no assigned routes");
183                this.A_er = netPlan.getMatrixLink2RouteAssignment();
184                this.A_se = DoubleFactory2D.dense.make(1 + nSRGs , Efm , 1.0); // 1 if link OK, 0 if fails
185                for (int contSRG = 0 ; contSRG < nSRGs ; contSRG ++)
186                        for (Link e : netPlan.getSRG (contSRG).getLinks())
187                                A_se.set(contSRG+1 , e.getIndex () , 0.0);
188                this.A_rs = DoubleFactory2D.dense.make(R , 1 + nSRGs , 1.0); // 1 if link OK, 0 if fails
189                for (Route r : netPlan.getRoutes ()) for (Link e : r.getSeqLinksRealPath()) for (SharedRiskGroup srg : e.getSRGs())
190                                        A_rs.set(r.getIndex (),1+srg.getIndex () , 0.0);
191
192                /* Check if the problem may have a solution: demands have at least one path in any failure state  */
193                if (A_dr.zMult(A_rs, null).getMinLocation() [0] == 0) { System.out.println("A_dr.zMult(A_rs, null): " + A_dr.zMult(A_rs, null)); throw new Net2PlanException ("Some demands cannot be routed in some failure state. We need more paths! (e.g. extend the path reach)"); }
194
195                this.Aout_ne = netPlan.getMatrixNodeLinkOutgoingIncidence();
196                this.Ain_ne = netPlan.getMatrixNodeLinkIncomingIncidence(); 
197                this.Abid_dd = DoubleFactory2D.sparse.make(D,D); for (Entry<Demand,Demand> entry : opposite_d.entrySet()) { Abid_dd.set(entry.getKey().getIndex (), entry.getValue().getIndex (), 1.0); Abid_dd.set(entry.getValue().getIndex (), entry.getKey().getIndex (), 1.0); } 
198                this.Abid_rr = DoubleFactory2D.sparse.make(R,R); for (Entry<Route,Route> entry : opposite_r.entrySet()) { Abid_rr.set(entry.getKey().getIndex (), entry.getValue().getIndex (), 1.0); Abid_rr.set(entry.getValue().getIndex (), entry.getKey().getIndex (), 1.0); } 
199                this.Abid_ee = DoubleFactory2D.sparse.make(Efm,Efm); for (Link e : netPlan.getLinks ()) { Abid_ee.set(e.getIndex (), oppositeLink(e).getIndex (), 1.0);  Abid_ee.set(oppositeLink(e).getIndex (), e.getIndex (), 1.0); }
200                
201                final long algorithmInitialtime = System.nanoTime();
202                final long algorithmEndtime = algorithmInitialtime + (long) (algorithm_maxExecutionTimeInSeconds.getDouble() * 1E9);
203
204                /* Up stage */
205                ArrayList<Integer> shuffledDemands = new ArrayList<Integer> (D); for (int d = 0 ; d < D ; d ++) shuffledDemands.add(d);
206                Collections.shuffle(shuffledDemands , rng);
207                if (tcfa_recoveryType.getString ().equals("shared"))
208                {
209                        DoubleMatrix1D best_xr = null; DoubleMatrix1D best_pe = null; double bestCost = Double.MAX_VALUE;
210                        int iterationCounter = 0;
211                        while ((System.nanoTime() < algorithmEndtime) && (iterationCounter < tcfa_maxNumIterations.getInt ()))
212                        {
213                                DoubleMatrix1D thisIt_xr = greedyAndLocalSearch_shared(algorithmEndtime , netPlan);
214                                DoubleMatrix1D thisIt_pe = A_er.zMult(thisIt_xr, null); thisIt_pe.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil);
215                                double thisItCost = computeCost_shared(thisIt_xr, thisIt_pe);
216                                if (thisItCost < bestCost) { best_xr = thisIt_xr; best_pe = thisIt_pe; bestCost = thisItCost; }
217                                iterationCounter ++;
218                        }
219                        stat_numItGRASP = iterationCounter;
220                        for (Route r : netPlan.getRoutes ()) r.setCarriedTraffic(tcfa_circuitCapacity_Gbps.getDouble () * best_xr.get(r.getIndex ()) , tcfa_circuitCapacity_Gbps.getDouble () * best_xr.get(r.getIndex ()));
221                        for (Link e : netPlan.getLinks ()) e.setCapacity(tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble () * best_pe.get(e.getIndex()));
222                } else if (tcfa_recoveryType.getString ().equals("1+1"))
223                {
224                        DoubleMatrix1D best_xr = null; DoubleMatrix1D best_x2r = null; DoubleMatrix1D best_pe = null; double bestCost = Double.MAX_VALUE;
225                        int iterationCounter = 0;
226                        while ((System.nanoTime() < algorithmEndtime) && (iterationCounter < tcfa_maxNumIterations.getInt ()))
227                        {
228                                Pair<DoubleMatrix1D,DoubleMatrix1D> thisItSolution = greedyAndLocalSearch_11(algorithmEndtime , netPlan);
229                                DoubleMatrix1D thisIt_xr = thisItSolution.getFirst();   
230                                DoubleMatrix1D thisIt_x2r = thisItSolution.getSecond();
231                                DoubleMatrix1D thisIt_pe = computePe_11(thisIt_xr, thisIt_x2r);
232                                double thisItCost = computeCost_11(thisIt_xr, thisIt_x2r , thisIt_pe);
233                                if (thisItCost < bestCost) { best_xr = thisIt_xr; best_x2r = thisIt_x2r; best_pe = thisIt_pe; bestCost = thisItCost; }
234                                iterationCounter ++;
235                        }
236                        stat_numItGRASP = iterationCounter;
237                        for (Demand d : netPlan.getDemands ())
238                        {
239                                if (Math.abs(best_xr.viewSelection(routeList_d [d.getIndex ()]).zSum() - 1) > 1E-3) throw new RuntimeException ("Bad");
240                                if (Math.abs(best_x2r.viewSelection(routeList_d [d.getIndex ()]).zSum() - 1) > 1E-3) throw new RuntimeException ("Bad");
241                                Route primaryRoute = null; ProtectionSegment protSegment = null;  
242                                for (Route r : d.getRoutes ()) if (Math.abs(best_xr.get(r.getIndex ()) - 1) <= 1e-3) { primaryRoute = r; primaryRoute.setCarriedTraffic(tcfa_circuitCapacity_Gbps.getDouble () , tcfa_circuitCapacity_Gbps.getDouble ()); break; }
243                                for (Route r : d.getRoutes ()) if (Math.abs(best_x2r.get(r.getIndex ()) - 1) <= 1e-3) { protSegment = netPlan.addProtectionSegment(r.getSeqLinksRealPath() , tcfa_circuitCapacity_Gbps.getDouble () , null); break; }
244                                primaryRoute.addProtectionSegment(protSegment);
245                        }
246                        for (Link e : netPlan.getLinks ()) e.setCapacity(tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble () * best_pe.get(e.getIndex()));
247                
248                } else if (tcfa_recoveryType.getString ().equals("restoration"))
249                {
250                        DoubleMatrix2D best_xrs = null; DoubleMatrix1D best_pe = null; double bestCost = Double.MAX_VALUE;
251                        int iterationCounter = 0;
252                        while ((System.nanoTime() < algorithmEndtime) && (iterationCounter < tcfa_maxNumIterations.getInt ()))
253                        {
254                                DoubleMatrix2D thisIt_xrs = greedyAndLocalSearch_restoration(algorithmEndtime , netPlan);
255                                DoubleMatrix1D thisIt_pe = computePe_restoration (thisIt_xrs); 
256                                double thisItCost = computeCost_shared(thisIt_xrs.viewColumn(0), thisIt_pe);
257                                if (thisItCost < bestCost) { best_xrs = thisIt_xrs; best_pe = thisIt_pe; bestCost = thisItCost; }
258                                iterationCounter ++;
259                        }
260                        stat_numItGRASP = iterationCounter;
261                        for (Route r : netPlan.getRoutes ()) r.setCarriedTraffic(tcfa_circuitCapacity_Gbps.getDouble () * best_xrs.get(r.getIndex (),0) , tcfa_circuitCapacity_Gbps.getDouble () * best_xrs.get(r.getIndex (),0));
262                        for (Link e : netPlan.getLinks ()) e.setCapacity(tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble () * best_pe.get(e.getIndex()));
263                } else throw new RuntimeException ("Bad");              
264                
265                /* Remove unused routes and links */
266                for (Route r : new HashSet<Route>(netPlan.getRoutes())) if (r.getCarriedTraffic() == 0) r.remove();
267                for (Link e : new HashSet<Link>(netPlan.getLinks())) if (e.getCapacity() == 0) e.remove();
268                
269                Quadruple<Double,Double,Double,Double> q = computeCost (netPlan);
270                stat_totalCost = q.getFirst();
271                stat_costinks = q.getSecond();
272                stat_costNodes = q.getThird();
273                stat_costCirc = q.getFourth();
274                stat_numLinks = (int) (netPlan.getVectorLinkCapacity().zSum () / (tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble ())); 
275                stat_numCirc = (int) Math.round(netPlan.getDemandTotalCarriedTraffic() / tcfa_circuitCapacity_Gbps.getDouble ()); if (tcfa_recoveryType.getString ().equals("1+1")) stat_numCirc *=2; 
276                for (Node n : netPlan.getNodes()) if (n.getOutgoingLinks().size() > 2) stat_numNodesDeg3 ++; else stat_numNodesDeg2 ++;
277                stat_totalTimeSecs = (System.nanoTime() - algorithmInitialtime)*1e-9;
278                checkSolution (netPlan);
279
280                double averageLinkUtilization = netPlan.getVectorLinkUtilizationIncludingProtectionSegments().zSum() / stat_numLinks;
281                final String fileNameStem = algorithm_outputFileNameRoot.getString() + "_" + tcfa_srgType.getString() + "_c" + tcfa_circuitCapacity_Gbps.getDouble () + "_t" + algorithm_maxExecutionTimeInSeconds.getDouble() + "_r" + tcfa_recoveryType.getString ();  
282                try 
283                {
284                        PrintWriter pw = new PrintWriter (new File (fileNameStem + "_allResults.txt"));
285                        pw.println(stat_totalCost );
286                        pw.println( stat_costinks );
287                        pw.println( stat_costNodes );
288                        pw.println( stat_costCirc );
289                        pw.println( stat_numLinks );
290                        pw.println( stat_numCirc );
291                        pw.println( stat_numNodesDeg2 );
292                        pw.println( stat_numNodesDeg3 );
293                        pw.println( stat_numItGRASP );
294                        pw.println( stat_numLSIterationsReducingCost );
295                        pw.println( stat_numSolverCalls );
296                        pw.println( stat_numSolverCallsOk );
297                        pw.println( stat_totalTimeSecs );
298                        pw.println( stat_totalTimeCreateProgramsSecs );
299                        pw.println( stat_totalTimeSolverSecs );
300                        pw.println( averageLinkUtilization );
301                        pw.close ();
302                } catch (Exception e) { e.printStackTrace(); throw new RuntimeException ("Not possible to write in File " + fileNameStem  + "_allResults.txt"); } 
303
304
305                netPlan.setAttribute("stat_totalCost", ""+stat_totalCost);
306                netPlan.setAttribute("stat_costinks", ""+stat_costinks);
307                netPlan.setAttribute("stat_costNodes", ""+stat_costNodes);
308                netPlan.setAttribute("stat_costCirc", ""+stat_costCirc);
309                netPlan.setAttribute("stat_numLinks", ""+stat_numLinks);
310                netPlan.setAttribute("stat_numCirc", ""+stat_numCirc);
311                netPlan.setAttribute("stat_numNodesDeg2", ""+stat_numNodesDeg2);
312                netPlan.setAttribute("stat_numNodesDeg3", ""+stat_numNodesDeg3);
313                netPlan.setAttribute("stat_numItGRASP", ""+stat_numItGRASP);
314                netPlan.setAttribute("stat_numLSIterationsReducingCost", ""+stat_numLSIterationsReducingCost);
315                netPlan.setAttribute("stat_numSolverCalls", ""+stat_numSolverCalls);
316                netPlan.setAttribute("stat_numSolverCallsOk", ""+stat_numSolverCallsOk);
317                netPlan.setAttribute("stat_totalTimeSecs", ""+stat_totalTimeSecs);
318                netPlan.setAttribute("stat_totalTimeCreateProgramsSecs", ""+stat_totalTimeCreateProgramsSecs);
319                netPlan.setAttribute("stat_totalTimeSolverSecs", ""+stat_totalTimeSolverSecs);
320
321                netPlan.saveToFile(new File (fileNameStem + ".n2p"));
322
323//              if (netPlan.getDemandTotalBlockedTraffic() > 1E-3) throw new RuntimeException ("Bad");
324//              if (netPlan.getLinksOversubscribed().size() > 0) throw new RuntimeException ("Bad");
325                System.out.println("JOM sol: cost: " + stat_totalCost + ", num links: "+ stat_numLinks);
326
327                return "Ok! cost: " + stat_totalCost + ", num links: "+ stat_numLinks;
328
329                } catch (Exception e) { e.printStackTrace(); throw new RuntimeException ("BAD OUT"); }  
330
331        }
332
333        @Override
334        public String getDescription()
335        {
336                return "This algorithm is devoted to solve the several network planning problems appearing in the case study in Section 12.10 of the book mentioned below. The context is the planning of a WDM optical network, given a set of node locations, and a set of end-to-end optical connections (lightpaths) to establish. The algorithm should decide on the fibers to deploy, and the routing of the lightpaths on the fibers, minimizing the total network cost. The cost model includes the fiber renting cost (paid to a dark fiber provider), OADM cost (separating between OADMs of degree two, and degree higher than two), and the cost of the transponders. The network must be tolerant to a set of user-defined failures (defined by SRGs). The algorithm permits choosing between two recovery types: 1+1, and shared protection. In the latter, the number of lightpaths between two nodes that survive to any of the failures is at least the number of optical connections requested for that node pair. For solving the problem, the algorithm uses a GRASP heuristic, to govern the iterative search that is performed using small ILPs solved with JOM in each iteration.";
337        }
338
339        @Override
340        public List<Triple<String, String, String>> getParameters()
341        {
342                /* Returns the parameter information for all the InputParameter objects defined in this object (uses Java reflection) */
343                return InputParameter.getInformationAllInputParameterFieldsOfObject(this);
344        }
345
346        private Quadruple<Double,Double,Double,Double> computeCost (NetPlan  np)
347        {
348                double costLinks = 0;
349                for (Link e : np.getLinks ())
350                        costLinks += linkCost_e.get(e.getIndex ()) * e.getCapacity() / (tcfa_linkCapacity_numCirc.getInt () * tcfa_circuitCapacity_Gbps.getDouble ());
351                
352                double costCircuits = 0;
353                for (Route r : np.getRoutes ())
354                        costCircuits += tcfa_circuitCost.getDouble () * (r.getCarriedTraffic() / tcfa_circuitCapacity_Gbps.getDouble ());
355                for (ProtectionSegment p : np.getProtectionSegments())
356                        costCircuits += tcfa_circuitCost.getDouble () * (p.getReservedCapacityForProtection() / tcfa_circuitCapacity_Gbps.getDouble ());
357                double costNodes = 0;
358                for (Node n : np.getNodes())
359                {
360                        final int nodeDegree = (int) Math.max(n.getOutgoingLinks().size(), n.getIncomingLinks().size()); 
361                        costNodes += (nodeDegree <= 2)? tcfa_costNodeDegree2.getDouble () : tcfa_costNodeDegreeHigherThan2.getDouble ();
362                }
363                final double totalCost = costLinks+  costCircuits + costNodes;
364
365                System.out.println("Total cost: " + totalCost + ", links %: " + (costLinks/totalCost) + ", circ %: " + (costCircuits/totalCost) + ", nodes %: " + (costNodes/totalCost));
366                
367                System.out.println("-- From np cost total :" + totalCost + ", links: " + costLinks + ", nodes: " + costNodes + ", circ: " + costCircuits);
368
369                return Quadruple.of(totalCost,costLinks,costNodes,costCircuits);
370        }
371        
372        private double computeCost_shared (DoubleMatrix1D x_r , DoubleMatrix1D p_e)
373        {
374                double costLinks = p_e.zDotProduct(linkCost_e);
375                final double costCircuits = tcfa_circuitCost.getDouble () * x_r.zSum();
376                double costNodes = 0;
377                DoubleMatrix1D degree_n = Aout_ne.zMult(p_e, null);
378                for (int n = 0; n < N ; n ++)
379                        costNodes += (degree_n.get(n) <= 2)? tcfa_costNodeDegree2.getDouble () : tcfa_costNodeDegreeHigherThan2.getDouble ();
380                final double totalCost = costLinks+  costCircuits + costNodes;
381                //System.out.println("-- shared cost total :" + totalCost + ", links: " + costLinks + ", nodes: " + costNodes + ", circ: " + costCircuits);
382                return totalCost;
383        }
384        private double computeCost_11 (DoubleMatrix1D x_r , DoubleMatrix1D x2_r , DoubleMatrix1D p_e)
385        {
386                double costLinks = p_e.zDotProduct(linkCost_e);
387                final double costCircuits = tcfa_circuitCost.getDouble () * (x_r.zSum() + x2_r.zSum());
388                double costNodes = 0;
389                DoubleMatrix1D degree_n = Aout_ne.zMult(p_e, null);
390                for (int n = 0; n < N ; n ++)
391                        costNodes += (degree_n.get(n) <= 2)? tcfa_costNodeDegree2.getDouble () : tcfa_costNodeDegreeHigherThan2.getDouble ();
392                final double totalCost = costLinks+  costCircuits + costNodes;
393                //System.out.println("-- shared cost total :" + totalCost + ", links: " + costLinks + ", nodes: " + costNodes + ", circ: " + costCircuits);
394                return totalCost;
395        }
396
397        private Link oppositeLink (Link e) { return e.getNetPlan ().getNodePairLinks(e.getDestinationNode(), e.getOriginNode() , false).iterator().next(); }
398        
399        private Pair<Map<Demand,Demand>,Map<Route,Route>> initializeNetPlanLinksBidirDemandsAndRoutes (NetPlan np)
400        {
401                /* Remove lower half demands from np */
402                np.removeAllRoutes(); np.removeAllProtectionSegments();
403                for (Node n1 : np.getNodes()) for (Node n2 : np.getNodes()) if (n1.getIndex () > n2.getIndex ()) for (Demand d : np.getNodePairDemands(n1, n2,false)) d.remove ();
404                np.addRoutesFromCandidatePathList(null , "K" , "" + tcfa_maxNumberPathsPerDemand.getInt () , "maxLengthInKm" , "" + tcfa_maxPathLengthInKm.getDouble() ,  "maxNumHops" , "" + tcfa_maxPathNumberOfHops.getInt ());
405
406                /* Add symmetric demands and routes */
407                Map<Demand,Demand> opposite_d = new HashMap<Demand,Demand> ();
408                Map<Route,Route> opposite_r = new HashMap<Route,Route> ();
409                for (Demand d : new HashSet<Demand> (np.getDemands()))
410                {
411                        final Demand opDemand = np.addDemand(d.getEgressNode(), d.getIngressNode(), d.getOfferedTraffic(), null);
412                        opposite_d.put(d,opDemand);
413                        opposite_d.put(opDemand,d);
414                        for (Route r : new HashSet<Route> (d.getRoutes ()))
415                        {
416                                final Route oppRoute = np.addRoute(opDemand, r.getCarriedTraffic(), r.getOccupiedCapacity() , oppositeSeqLinks (r.getSeqLinksRealPath()), null);
417                                opposite_r.put(r,oppRoute); opposite_r.put(oppRoute,r);
418                        }
419                }
420                return Pair.of(opposite_d,opposite_r);
421        }
422
423        private void initializeSRGs (NetPlan np)
424        {
425                np.removeAllSRGs();
426                if (tcfa_srgType.getString ().equals("perBidirectionalLinkBundle"))
427                {
428                        for (Node n1 : np.getNodes()) for (Node n2 : np.getNodes()) if (n1.getIndex () < n2.getIndex())
429                        {
430                                final double linkLengthKm = np.getNodePairLinks(n1, n2,false).iterator().next().getLengthInKm();
431                                final SharedRiskGroup srg = np.addSRG(tcfa_srgMttfPer1000Km_hrs.getDouble () * linkLengthKm / 1000, tcfa_srgMttr_hrs.getDouble (), null);
432                                for (Link e : np.getNodePairLinks(n1, n2,true)) srg.addLink(e); 
433                        }
434                }
435                else if (tcfa_srgType.getString ().equals("noFailure"))
436                {
437                } else throw new Net2PlanException ("Wrong SRG type");
438        }
439        
440        private List<Link> oppositeSeqLinks  (List<Link> s)
441        {
442                LinkedList<Link> oppList = new LinkedList<Link> ();
443                for (Link e : s) oppList.addFirst(oppositeLink(e));
444                return (List<Link>) oppList;
445        }
446
447        private void checkSolution (NetPlan np)
448        {
449                if (np.getDemandsBlocked().size() != 0)
450                {
451                        System.out.println("np.getDemandsBlocked(): " + np.getDemandsBlocked());
452                        System.out.println("np.getDemandOfferedTrafficMap(): " + np.getVectorDemandOfferedTraffic());
453                        System.out.println("np.getDemandBlockedTrafficMap(): " + np.getVectorDemandBlockedTraffic());
454                        throw new RuntimeException ("Bad");
455                }
456                if (np.getLinksOversubscribed().size() != 0)
457                {
458                        System.out.println("np.getLinksOversubscribed(): " + np.getLinksOversubscribed());
459                        System.out.println("np.getLinkCapacityMap(): " + np.getVectorLinkCapacity());
460                        System.out.println("np.getLinkCarriedTrafficMap(): " + np.getVectorLinkTotalCarriedTraffic());
461                        throw new RuntimeException ("Bad");
462                }
463                
464                /* Check bidirectional topology */
465                for (Node n1 : np.getNodes()) 
466                        for (Node n2 : np.getNodes()) 
467                                if (n1.getIndex () > n2.getIndex ()) 
468                                {
469                                        Set<Link> thisLinks = np.getNodePairLinks(n1, n2,false);
470                                        Set<Link> oppLinks = np.getNodePairLinks(n2, n1,false);
471                                        if (thisLinks.size() > 1) throw new RuntimeException ("Bad");
472                                        if (oppLinks.size() > 1) throw new RuntimeException ("Bad");
473                                        if (tcfa_bidirectionalLinks.getBoolean())
474                                        {
475                                                if (thisLinks.size() != oppLinks.size()) throw new RuntimeException ("Bad");
476                                                if (thisLinks.size() > 0)
477                                                        if (thisLinks.iterator().next().getCapacity() != oppLinks.iterator().next().getCapacity() ) throw new RuntimeException ("Bad");
478                                        }
479                                }
480
481                /* Check bidirectional routes */
482                for (Route r : np.getRoutes())
483                {
484                        final Route oppRoute = opposite_r.get(r);
485                        if (!np.getRoutes().contains(oppRoute)) throw new RuntimeException ("Bad");
486                        if (r.getCarriedTraffic() != oppRoute.getCarriedTraffic()) throw new RuntimeException ("Bad");
487                }
488
489                if (tcfa_recoveryType.getString ().equals("1+1"))
490                { // one prot segment, and link disjoint
491                        for (Route r : np.getRoutes())
492                        {
493                                if (r.getPotentialBackupProtectionSegments().size() != 1) throw new RuntimeException ("Bad");
494                                final ProtectionSegment protSegment = r.getPotentialBackupProtectionSegments().iterator().next();
495                                List<Link> seqLinks = new ArrayList<Link> (r.getSeqLinksRealPath());
496                                seqLinks.retainAll(protSegment.getSeqLinks());
497                                if (!seqLinks.isEmpty()) 
498                                {
499                                        System.out.println ("Route : "+  r + " links: " + r.getSeqLinksRealPath() + ", prot segment " + protSegment + " links: " + protSegment.getSeqLinks() + ", get route carried traffic: " + r.getCarriedTraffic() + ", prot segment reserved traffic: " + protSegment.getReservedCapacityForProtection());
500                                        throw new RuntimeException ("Bad");
501                                }
502                        }
503                }
504                
505        }
506        
507        private DoubleMatrix1D jomStep_shared (Set<Integer> demandsChange , DoubleMatrix1D maxPe , DoubleMatrix1D  x_r)
508        {
509                final long initTimeBeforeCreatingProgram = System.nanoTime();
510                if (demandsChange.isEmpty()) throw new RuntimeException ("Bad");
511                final int [] dChange = IntUtils.toArray(demandsChange);
512                int [] rChange = new int [0]; for (int d : dChange) rChange = IntUtils.concatenate(rChange,  routeList_d [d]);
513                final int Rrest = rChange.length;
514                x_r.viewSelection(rChange).assign(0);
515                DoubleMatrix2D A_dRestrRest = A_dr.viewSelection(dChange, rChange);
516                DoubleMatrix2D A_erRest = A_er.viewSelection(null, rChange);
517                DoubleMatrix2D A_rRests = A_rs.viewSelection(rChange,null);
518                DoubleMatrix1D occupiedWithoutChangingR_e = A_er.zMult (x_r,null);
519                DoubleMatrix2D Abid_rRestrRest = DoubleFactory2D.sparse.make(Rrest,Rrest);
520                for (int rRest = 0 ; rRest < Rrest ; rRest ++) 
521                { 
522                        final int r = rChange [rRest]; 
523                        final int opp_r = opposite_r.get(netPlan.getRoute(r)).getIndex ();
524                        final int [] oppRrest = IntUtils.find(rChange, opp_r, SearchType.ALL);
525                        if (oppRrest.length != 1) throw new RuntimeException ("Bad");
526                        Abid_rRestrRest.set(rRest , oppRrest [0] , 1.0); Abid_rRestrRest.set(oppRrest [0], rRest , 1);
527                }
528                OptimizationProblem op = new OptimizationProblem();
529                op.setInputParameter("R",R);
530                op.setInputParameter("U",tcfa_linkCapacity_numCirc.getInt ());
531                op.setInputParameter("c_e", linkCost_e.toArray() , "row");
532                op.setInputParameter("c",  tcfa_circuitCost.getDouble ());
533                op.setInputParameter("d3Cost",  tcfa_costNodeDegreeHigherThan2.getDouble ());
534                op.setInputParameter("onesS",  DoubleUtils.ones(1+nSRGs) , "row");
535                op.setInputParameter("onesE",  DoubleUtils.ones(Efm) , "row");
536                op.setInputParameter("A_dRestrRest", new DoubleMatrixND (A_dRestrRest));
537                op.setInputParameter("h_dRest", numCircH_d.viewSelection(dChange).toArray(), "row");
538                op.setInputParameter("A_rRests", new DoubleMatrixND (A_rRests));
539                op.setInputParameter("A_erRest", new DoubleMatrixND (A_erRest));
540                op.setInputParameter("Aout_ne", new DoubleMatrixND (Aout_ne));
541                op.setInputParameter("Abid_rRestrRest", new DoubleMatrixND (Abid_rRestrRest));
542                op.setInputParameter("occup_e", occupiedWithoutChangingR_e.toArray(),  "row");
543                
544                op.addDecisionVariable("d3_n", true, new int[] { 1, N }, 0, 1); 
545                op.addDecisionVariable("p_e", true, new int[] { 1, Efm }, DoubleUtils.zeros(Efm) , (maxPe == null)? DoubleUtils.constantArray(Efm, Double.MAX_VALUE) : maxPe.toArray());
546                op.addDecisionVariable("x_rRest", true, new int[] { 1, Rrest}, 0 ,  Double.MAX_VALUE); /* 1 if there is a link from node i to node j, 0 otherwise */
547
548                op.setObjectiveFunction("minimize", "sum(c_e .* p_e) + c * sum (x_rRest) + d3Cost * sum(d3_n)");
549                op.addConstraint("A_dRestrRest * diag(x_rRest) * A_rRests >= h_dRest' * onesS"); /* the flow-conservation constraints (NxD constraints) */
550                op.addConstraint("A_erRest * x_rRest' + occup_e' <= U * p_e'"); /* the capacity constraints (E constraints) */
551                op.addConstraint("Aout_ne * p_e' <= 2 + R * d3_n'"); /* the capacity constraints (E constraints) */
552                op.addConstraint("x_rRest == x_rRest * Abid_rRestrRest");
553                if (tcfa_bidirectionalLinks.getBoolean()) { op.setInputParameter("Abid_ee", new DoubleMatrixND (Abid_ee)); op.addConstraint("p_e == p_e * Abid_ee"); } 
554
555                stat_totalTimeCreateProgramsSecs += eTime (initTimeBeforeCreatingProgram); 
556                                
557                double solverTimeAllowed = algorithm_maxSolverTimeInSeconds.getDouble();
558                do 
559                {
560                        final long t = System.nanoTime(); 
561                        try 
562                        { 
563                                stat_numSolverCalls ++; op.solve(algorithm_solverName.getString() , "solverLibraryName", algorithm_solverLibraryName.getString () ,  "maxSolverTimeInSeconds" , solverTimeAllowed);
564                        } catch (Exception e) { System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); solverTimeAllowed += 2; stat_totalTimeSolverSecs += eTime(t);  continue;} 
565                        stat_totalTimeSolverSecs += eTime(t);  
566                        if (solverTimeAllowed != algorithm_maxSolverTimeInSeconds.getDouble()) System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed);
567                        solverTimeAllowed *= 2; // duplicate the time if not feasible solution is found
568                } while (!op.solutionIsFeasible());             
569                stat_numSolverCallsOk ++;
570                final double [] new_x_rRest = op.getPrimalSolution("x_rRest").to1DArray();
571                x_r.viewSelection(rChange).assign(new_x_rRest);
572                return x_r;
573        }
574
575        private Pair<DoubleMatrix1D,DoubleMatrix1D> jomStep_11 (Set<Integer> demandsChange ,  DoubleMatrix1D maxPe , DoubleMatrix1D x_r , DoubleMatrix1D x2_r)
576        {
577//              System.out.println ("demandsChange :" + demandsChange);
578                final long initTimeBeforeCreatingProgram = System.nanoTime();
579
580                if (demandsChange.isEmpty()) throw new RuntimeException ("Bad");
581                final int [] dChange = IntUtils.toArray(demandsChange);
582//              System.out.println ("dChange :" + Arrays.toString(dChange));
583                int [] rChange = new int [0]; for (int d : dChange) rChange = IntUtils.concatenate(rChange,  routeList_d [d]);
584                final int Rrest = rChange.length;
585                final int Drest = dChange.length;
586                x_r.viewSelection(rChange).assign(0);
587                x2_r.viewSelection(rChange).assign(0);
588                DoubleMatrix2D A_dRestrRest = A_dr.viewSelection(dChange, rChange);
589                DoubleMatrix2D A_erRest = A_er.viewSelection(null, rChange);
590                DoubleMatrix2D A_rRests = A_rs.viewSelection(rChange,null);
591                DoubleMatrix1D occupiedWithoutChangingR_e = A_er.zMult (x_r,null);
592                DoubleMatrix1D backup_occupiedWithoutChangingR_e = A_er.zMult (x2_r,null);
593                for (int e = 0;e < Efm; e ++) occupiedWithoutChangingR_e.set(e, occupiedWithoutChangingR_e.get(e) + backup_occupiedWithoutChangingR_e.get(e));
594                DoubleMatrix2D Abid_rRestrRest = DoubleFactory2D.sparse.make(Rrest,Rrest);
595                for (int rRest = 0 ; rRest < Rrest ; rRest ++) 
596                { 
597                        final int r = rChange [rRest]; final int opp_r = opposite_r.get(netPlan.getRoute (r)).getIndex ();
598                        final int [] oppRrest = IntUtils.find(rChange, opp_r, SearchType.ALL); if (oppRrest.length != 1) throw new RuntimeException ("Bad");
599                        Abid_rRestrRest.set(rRest , oppRrest [0] , 1.0); Abid_rRestrRest.set(oppRrest [0], rRest , 1);
600                }
601                
602                OptimizationProblem op = new OptimizationProblem();
603                op.setInputParameter("R",R);
604                op.setInputParameter("U",tcfa_linkCapacity_numCirc.getInt ());
605                op.setInputParameter("c_e", linkCost_e.toArray() , "row");
606                op.setInputParameter("c",  tcfa_circuitCost.getDouble ());
607                op.setInputParameter("d3Cost",  tcfa_costNodeDegreeHigherThan2.getDouble ());
608                op.setInputParameter("onesS",  DoubleUtils.ones(1+nSRGs) , "row");
609                op.setInputParameter("onesE",  DoubleUtils.ones(Efm) , "row");
610                op.setInputParameter("A_dRestrRest", new DoubleMatrixND (A_dRestrRest));
611                op.setInputParameter("h_dRest", numCircH_d.viewSelection(dChange).toArray(), "row");
612                op.setInputParameter("A_rRests", new DoubleMatrixND (A_rRests));
613                op.setInputParameter("A_erRest", new DoubleMatrixND (A_erRest));
614                op.setInputParameter("Aout_ne", new DoubleMatrixND (Aout_ne));
615                op.setInputParameter("Abid_rRestrRest", new DoubleMatrixND (Abid_rRestrRest));
616                op.setInputParameter("occup_e", occupiedWithoutChangingR_e.toArray(),  "row");
617                
618                op.addDecisionVariable("d3_n", true, new int[] { 1, N }, 0, 1); 
619                op.addDecisionVariable("p_e", true, new int[] { 1, Efm }, DoubleUtils.zeros(Efm) , (maxPe == null)? DoubleUtils.constantArray(Efm, Double.MAX_VALUE) : maxPe.toArray());
620                op.addDecisionVariable("x_rRest", true, new int[] { 1, Rrest}, 0 ,  1); 
621                op.addDecisionVariable("x2_rRest", true, new int[] { 1, Rrest}, 0 ,  1); 
622
623                op.setObjectiveFunction("minimize", "sum(c_e .* p_e) + d3Cost * sum(d3_n)");
624                op.addConstraint("A_dRestrRest * x_rRest' == h_dRest'"); /* the flow-conservation constraints (NxD constraints) */
625                op.addConstraint("A_dRestrRest * x2_rRest' == h_dRest'"); /* the flow-conservation constraints (NxD constraints) */
626                op.addConstraint("A_erRest * (x_rRest+x2_rRest)' + occup_e' <= U * p_e'"); /* the capacity constraints (E constraints) */
627                op.addConstraint("Aout_ne * p_e' <= 2 + R*d3_n'"); /* the capacity constraints (E constraints) */
628                op.addConstraint("A_dRestrRest * diag(x_rRest+x2_rRest) * A_erRest' <= 1"); 
629                op.addConstraint("x_rRest == x_rRest * Abid_rRestrRest"); 
630                op.addConstraint("x2_rRest == x2_rRest * Abid_rRestrRest"); 
631                if (tcfa_bidirectionalLinks.getBoolean()) { op.setInputParameter("Abid_ee", new DoubleMatrixND (Abid_ee)); op.addConstraint("p_e == p_e * Abid_ee"); } 
632
633//              System.out.println("Abid_rr sum filas: " + Abid_rRestrRest.zMult(DoubleFactory1D.dense.make(Rrest,1), null));
634//              System.out.println("Abid_rr sum cols: " +  Abid_rRestrRest.zMult(DoubleFactory1D.dense.make(Rrest,1), null , 1 , 0 , true));
635//              System.out.println("Abid_rr non zeros: " +  Abid_rRestrRest.zSum());
636                stat_totalTimeCreateProgramsSecs += eTime (initTimeBeforeCreatingProgram); 
637                
638                double solverTimeAllowed = algorithm_maxSolverTimeInSeconds.getDouble ();
639                do 
640                { 
641                        final long t = System.nanoTime(); 
642                        try 
643                        {
644                                stat_numSolverCalls ++; op.solve(algorithm_solverName.getString() , "solverLibraryName", algorithm_solverLibraryName.getString () ,  "maxSolverTimeInSeconds" , solverTimeAllowed);
645                        } catch (Exception e) { System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); solverTimeAllowed += 2; stat_totalTimeSolverSecs += eTime(t);  continue; }
646                        if (solverTimeAllowed != algorithm_maxSolverTimeInSeconds.getDouble ()) System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed);
647                        stat_totalTimeSolverSecs += eTime(t);  
648                        solverTimeAllowed *= 2; // duplicate the time if not feasible solution is found
649                        //System.out.println ("Time to create program: " + (1E-9*(t-initTimeBeforeCreatingProgram)) + ", solver time " + ((System.nanoTime()-t)*1E-9) + " s ; ");
650                } while (!op.solutionIsFeasible());             
651                stat_numSolverCallsOk ++; 
652                final double[] new_x_rRest = op.getPrimalSolution("x_rRest").to1DArray();
653                final double[] new_x2_rRest = op.getPrimalSolution("x2_rRest").to1DArray();
654
655                x_r.viewSelection(rChange).assign(new_x_rRest);
656                x2_r.viewSelection(rChange).assign(new_x2_rRest);
657                return Pair.of(x_r,x2_r);
658        }
659
660        private DoubleMatrix2D jomStep_restoration (Set<Integer> demandsChange ,  DoubleMatrix1D maxPe , DoubleMatrix2D  x_rs)
661        {
662                final long initTimeBeforeCreatingProgram = System.nanoTime();
663                if (demandsChange.isEmpty()) throw new RuntimeException ("Bad");
664                final int [] dChange = IntUtils.toArray(demandsChange);
665                int [] rChange = new int [0]; for (int d : dChange) rChange = IntUtils.concatenate(rChange,  routeList_d [d]);
666                final int Rrest = rChange.length;
667                final int Drest = dChange.length;
668                x_rs.viewSelection(rChange,null).assign(0);
669//              DoubleMatrix2D A_dRestrRest = A_dr.viewSelection(dChange, rChange);
670                DoubleMatrix2D A_dRestrRest = DoubleFactory2D.sparse.make(Drest,Rrest);
671                for (int rRest = 0 ; rRest < Rrest ; rRest ++) 
672                { 
673                        final int r = rChange [rRest]; final int dRest [] = IntUtils.find(dChange, netPlan.getRoute (r).getDemand().getIndex () , SearchType.ALL);
674                        if (dRest.length  != 1) throw new RuntimeException ("Bad");
675                        A_dRestrRest.set(dRest [0], rRest , 1.0);
676                }
677                
678                DoubleMatrix2D A_erRest = A_er.viewSelection(null, rChange);
679                DoubleMatrix2D A_rRests = A_rs.viewSelection(rChange,null);
680                DoubleMatrix2D occupiedWithoutChangingR_es = A_er.zMult (x_rs,null);
681                DoubleMatrix2D Abid_rRestrRest = DoubleFactory2D.sparse.make(Rrest,Rrest);
682                for (int rRest = 0 ; rRest < Rrest ; rRest ++) 
683                { 
684                        final int r = rChange [rRest]; final int opp_r = opposite_r.get(netPlan.getRoute (r)).getIndex ();
685                        final int [] oppRrest = IntUtils.find(rChange, opp_r, SearchType.ALL);
686                        if (oppRrest.length != 1) throw new RuntimeException ("Bad");
687                        Abid_rRestrRest.set(rRest , oppRrest [0] , 1.0); Abid_rRestrRest.set(oppRrest [0], rRest , 1);
688                }
689                
690                OptimizationProblem op = new OptimizationProblem();
691                op.setInputParameter("R",R);
692                op.setInputParameter("U",tcfa_linkCapacity_numCirc.getInt ());
693                op.setInputParameter("c_e", linkCost_e.toArray() , "row");
694                op.setInputParameter("c",  tcfa_circuitCost.getDouble ());
695                op.setInputParameter("d3Cost",  tcfa_costNodeDegreeHigherThan2.getDouble ());
696                op.setInputParameter("onesS",  DoubleUtils.ones(1+nSRGs) , "row");
697                op.setInputParameter("onesE",  DoubleUtils.ones(Efm) , "row");
698                op.setInputParameter("A_dRestrRest", new DoubleMatrixND (A_dRestrRest));
699                op.setInputParameter("h_dRest", numCircH_d.viewSelection(dChange).toArray(), "row");
700                op.setInputParameter("A_rRests", new DoubleMatrixND (A_rRests));
701                op.setInputParameter("A_erRest", new DoubleMatrixND (A_erRest));
702                op.setInputParameter("Aout_ne", new DoubleMatrixND (Aout_ne));
703                op.setInputParameter("Abid_rRestrRest", new DoubleMatrixND (Abid_rRestrRest));
704                op.setInputParameter("occup_es", occupiedWithoutChangingR_es.toArray());
705                
706                op.addDecisionVariable("d3_n", true, new int[] { 1, N }, 0, 1); 
707                op.addDecisionVariable("p_e", true, new int[] { 1, Efm }, DoubleUtils.zeros(Efm) , (maxPe == null)? DoubleUtils.constantArray(Efm, Double.MAX_VALUE) : maxPe.toArray());
708                op.addDecisionVariable("x_rRests", true, new int[] { Rrest, 1+nSRGs}, 0 , Double.MAX_VALUE); /* 1 if there is a link from node i to node j, 0 otherwise */
709
710                op.setObjectiveFunction("minimize", "sum(c_e .* p_e) + d3Cost * sum(d3_n)");
711                op.addConstraint("A_dRestrRest * x_rRests(all,0) == h_dRest'"); /* the flow-conservation constraints (NxD constraints) */
712                op.addConstraint("A_dRestrRest * (x_rRests .* A_rRests) >= h_dRest' * onesS"); /* the flow-conservation constraints (NxD constraints) */
713                op.addConstraint("A_erRest * x_rRests + occup_es <= U * p_e' * onesS"); /* the capacity constraints (E constraints) */
714                op.addConstraint("Aout_ne * p_e' <= 2 + R*d3_n'"); /* the capacity constraints (E constraints) */
715                op.addConstraint("x_rRests == Abid_rRestrRest * x_rRests"); 
716                op.addConstraint("x_rRests >= diag(x_rRests(all,0)) * A_rRests"); // if not affected by failure, do not change 
717                if (tcfa_bidirectionalLinks.getBoolean()) { op.setInputParameter("Abid_ee", new DoubleMatrixND (Abid_ee)); op.addConstraint("p_e == p_e * Abid_ee"); } 
718
719                stat_totalTimeCreateProgramsSecs += eTime (initTimeBeforeCreatingProgram); 
720
721                double solverTimeAllowed = algorithm_maxSolverTimeInSeconds.getDouble ();
722                do 
723                { 
724                        final long t = System.nanoTime(); 
725                        try { stat_numSolverCalls ++; op.solve(algorithm_solverName.getString() , "solverLibraryName", algorithm_solverLibraryName.getString () ,  "maxSolverTimeInSeconds" , solverTimeAllowed); } catch (Exception e) { System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed); solverTimeAllowed += 2; stat_totalTimeSolverSecs += eTime(t);  continue;} 
726                        if (solverTimeAllowed != algorithm_maxSolverTimeInSeconds.getDouble ()) System.out.println ("-- EXTRA TIME FOR SOLVER, NOW: " + solverTimeAllowed);
727                        stat_totalTimeSolverSecs += eTime(t);  
728                        solverTimeAllowed *= 2; // duplicate the time if not feasible solution is found
729                        //System.out.println ("Time to create program: " + (1E-9*(t-initTimeBeforeCreatingProgram)) + ", solver time " + ((System.nanoTime()-t)*1E-9) + " s ; ");
730                } while (!op.solutionIsFeasible());             
731                stat_numSolverCallsOk ++; 
732                final DoubleMatrix2D new_x_rRests = op.getPrimalSolution("x_rRests").view2D();
733                final DoubleMatrix1D new_pe = DoubleFactory1D.dense.make(op.getPrimalSolution("p_e").to1DArray());
734                x_rs.viewSelection(rChange, null).assign(new_x_rRests);
735                return x_rs;
736        }
737
738        private DoubleMatrix1D greedyAndLocalSearch_shared (long algorithmEndtime , NetPlan np)
739        {
740                final long initGreedy = System.nanoTime();
741                ArrayList<Pair<Node,Node>> shuffledNodePairs = new ArrayList<Pair<Node,Node>> (N*(N-1)/2); 
742                for (Node n1 : np.getNodes ()) for (Node n2 : np.getNodes ()) if (n1.getIndex () < n2.getIndex ()) shuffledNodePairs.add(Pair.of(n1, n2));
743                Collections.shuffle(shuffledNodePairs , rng);
744                DoubleMatrix1D current_xr = DoubleFactory1D.dense.make(R);
745                for (Pair<Node,Node> nodePair : shuffledNodePairs)
746                {
747                        Set<Integer> thisDemand = new HashSet<Integer> ();
748                        for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true))
749                                thisDemand.add(d.getIndex());
750//                      System.out.println("thsiDemand: " + thisDemand);
751                        if (jomStep_shared(thisDemand , null , current_xr) == null) throw new RuntimeException ("Returned null");
752                }
753                DoubleMatrix1D current_pe = A_er.zMult(current_xr, null); current_pe.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil);
754                double current_cost = computeCost_shared(current_xr, current_pe); 
755                //DoubleMatrix1D best_xr = current_xr.copy(); DoubleMatrix1D best_pe = current_pe.copy(); double bestCost = ;
756                
757                System.out.println("---- Greedy solution Time " + (1E-9*(System.nanoTime() - initGreedy)) +",  Cost : "+  current_cost + ", num links: " + current_pe.zSum());
758
759                while (System.nanoTime() < algorithmEndtime) 
760                {
761                        //DoubleMatrix1D bestNeighbor_xr = null; DoubleMatrix1D bestNeighbor_pe = null; double bestNeighborCost = Double.MAX_VALUE; Pair<Long,Long> bestNeighborNodePair = null;
762                        Collections.shuffle(shuffledNodePairs ,rng);
763                        final long initLSIteration = System.nanoTime();
764
765                        boolean improvedSolution = false;
766                        for (Pair<Node,Node> nodePair : shuffledNodePairs)
767                        {
768//                              System.out.println("Pair: " + nodePair + ". Tabu list: " + tabuList_n1n2 + ", IS TABU: " + tabuList_n1n2.contains(nodePair));
769                                Set<Integer> demandsToChange = new HashSet<Integer> (); 
770                                for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true))
771                                        demandsToChange.add(d.getIndex ()); 
772                                DoubleMatrix1D thisNeighbor_xr = current_xr.copy(); // do not modify the original one
773                                if (jomStep_shared(demandsToChange , null , thisNeighbor_xr) == null) throw new RuntimeException ("Returned null");
774                                DoubleMatrix1D thisNeighbor_pe = A_er.zMult(thisNeighbor_xr, null); thisNeighbor_pe.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil);
775                                final double thisNeighborCost = computeCost_shared (thisNeighbor_xr , thisNeighbor_pe);
776                                if (thisNeighborCost < current_cost)
777                                { 
778                                        current_xr = thisNeighbor_xr; current_pe = thisNeighbor_pe; current_cost = thisNeighborCost; improvedSolution = true; 
779                                        System.out.println("NEW BEST COST: " + current_cost + ", num links: " + current_pe.zSum()); 
780                                        break; 
781                                }
782                        }                       
783                        
784                        System.out.println("- Local search TIME: " + ((System.nanoTime()-initLSIteration)*1e-9) + ", current cost: " + current_cost + ", num links: " + current_pe.zSum() + ", IMPROVED: " + improvedSolution);
785
786                        if (!improvedSolution) break; else stat_numLSIterationsReducingCost ++; // end grasp iteration 
787                }
788                
789                return current_xr;
790        }
791
792        private Pair<DoubleMatrix1D,DoubleMatrix1D> greedyAndLocalSearch_11 (long algorithmEndtime , NetPlan np)
793        {
794                final long initGreedy = System.nanoTime();
795                ArrayList<Pair<Node,Node>> shuffledNodePairs = new ArrayList<Pair<Node,Node>> (N*(N-1)/2); 
796                for (Node n1 : np.getNodes ()) for (Node n2 : np.getNodes ()) if (n1.getIndex () < n2.getIndex ()) shuffledNodePairs.add(Pair.of(n1, n2));
797                Collections.shuffle(shuffledNodePairs , rng);
798                DoubleMatrix1D current_xr = DoubleFactory1D.dense.make(R);
799                DoubleMatrix1D current_x2r = DoubleFactory1D.dense.make(R);
800                for (Pair<Node,Node> nodePair : shuffledNodePairs)
801                {
802                        Set<Integer> thisDemand = new HashSet<Integer> ();
803                        for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true))
804                                thisDemand.add(d.getIndex());
805                        final long t = System.nanoTime(); if (jomStep_11(thisDemand , null , current_xr , current_x2r) == null) throw new RuntimeException ("Returned null");
806                        //System.out.println("Node pair " + nodePair + " time increasing it: " + (1E-9*(System.nanoTime()-t)) + ", num demands: " + thisDemand.size());
807                }
808                DoubleMatrix1D current_pe = computePe_11 (current_xr , current_xr);
809                double current_cost = computeCost_11(current_xr, current_x2r , current_pe); 
810
811                System.out.println("---- Greedy solution Time " + (1E-9*(System.nanoTime() - initGreedy)) +",  Cost : "+  current_cost + ", num links: " + current_pe.zSum());
812
813                while (System.nanoTime() < algorithmEndtime) 
814                {
815                        Collections.shuffle(shuffledNodePairs ,rng);
816                        final long initLSIteration = System.nanoTime();
817
818                        boolean improvedSolution = false;
819                        for (Pair<Node,Node> nodePair : shuffledNodePairs)
820                        {
821                                Set<Integer> demandsToChange = new HashSet<Integer> (); 
822                                for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true))
823                                        demandsToChange.add(d.getIndex()); 
824                                DoubleMatrix1D thisNeighbor_xr = current_xr.copy(); // do not modify the original one
825                                DoubleMatrix1D thisNeighbor_x2r = current_x2r.copy(); // do not modify the original one
826                                if (jomStep_11(demandsToChange , null , thisNeighbor_xr , thisNeighbor_x2r) == null) throw new RuntimeException ("Returned null");
827                                DoubleMatrix1D thisNeighbor_pe = computePe_11(thisNeighbor_xr, thisNeighbor_x2r);
828                                final double thisNeighborCost = computeCost_11 (thisNeighbor_xr , thisNeighbor_x2r , thisNeighbor_pe);
829                                if (thisNeighborCost < current_cost) 
830                                { 
831                                        current_xr = thisNeighbor_xr; current_x2r = thisNeighbor_x2r; current_pe = thisNeighbor_pe; current_cost = thisNeighborCost; improvedSolution = true; 
832                                        System.out.println("NEW BEST COST: " + current_cost + ", num links: " + current_pe.zSum()); 
833                                        break;
834                                }
835                        }                       
836                        
837                        System.out.println("- Local search TIME: " + ((System.nanoTime()-initLSIteration)*1e-9) + ", current cost: " + current_cost + ", num links: " + current_pe.zSum() + ", IMPROVED: " + improvedSolution);
838
839                        if (!improvedSolution) break; else stat_numLSIterationsReducingCost ++; // end grasp iteration 
840                }
841                
842                return Pair.of(current_xr , current_x2r);
843        }
844        private DoubleMatrix2D greedyAndLocalSearch_restoration (long algorithmEndtime , NetPlan np)
845        {
846                final long initGreedy = System.nanoTime();
847                ArrayList<Pair<Node,Node>> shuffledNodePairs = new ArrayList<Pair<Node,Node>> (N*(N-1)/2); 
848                for (Node n1 : np.getNodes ()) for (Node n2 : np.getNodes ()) if (n1.getIndex () < n2.getIndex()) shuffledNodePairs.add(Pair.of(n1, n2));
849                Collections.shuffle(shuffledNodePairs , rng);
850                DoubleMatrix2D current_xrs = DoubleFactory2D.dense.make(R,1+nSRGs);
851                for (Pair<Node,Node> nodePair : shuffledNodePairs)
852                {
853                        Set<Integer> thisDemand = new HashSet<Integer> ();
854                        for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond(),true))
855                                thisDemand.add(d.getIndex ());
856//                      System.out.println("thsiDemand: " + thisDemand);
857                        if (jomStep_restoration(thisDemand , null , current_xrs) == null) throw new RuntimeException ("Returned null");
858                }
859                DoubleMatrix1D current_pe = computePe_restoration(current_xrs);
860                double current_cost = computeCost_shared(current_xrs.viewColumn(0), current_pe); 
861                
862                System.out.println("---- Greedy solution Time " + (1E-9*(System.nanoTime() - initGreedy)) +",  Cost : "+  current_cost + ", num links: " + current_pe.zSum() + ", SRGs: " + nSRGs);
863
864                while (System.nanoTime() < algorithmEndtime) 
865                {
866                        Collections.shuffle(shuffledNodePairs ,rng);
867                        final long initLSIteration = System.nanoTime();
868
869                        boolean improvedSolution = false;
870                        for (Pair<Node,Node> nodePair : shuffledNodePairs)
871                        {
872                                Set<Integer> demandsToChange = new HashSet<Integer> (); 
873                                for (Demand d : np.getNodePairDemands(nodePair.getFirst(), nodePair.getSecond() , true))
874                                        demandsToChange.add(d.getIndex()); 
875                                DoubleMatrix2D thisNeighbor_xrs = current_xrs.copy(); // do not modify the original one
876                                if (jomStep_restoration(demandsToChange , null , thisNeighbor_xrs) == null) throw new RuntimeException ("Returned null");
877                                DoubleMatrix1D thisNeighbor_pe = computePe_restoration(thisNeighbor_xrs);
878                                final double thisNeighborCost = computeCost_shared (thisNeighbor_xrs.viewColumn(0) , thisNeighbor_pe);
879                                if (thisNeighborCost < current_cost) 
880                                { 
881                                        current_xrs = thisNeighbor_xrs; current_pe = thisNeighbor_pe; current_cost = thisNeighborCost; improvedSolution = true; 
882                                        System.out.println("NEW BEST COST: " + current_cost + ", num links: " + current_pe.zSum()); 
883                                        break; 
884                                }
885                        }                       
886                        
887                        System.out.println("- Local search TIME: " + ((System.nanoTime()-initLSIteration)*1e-9) + ", current cost: " + current_cost + ", num links: " + current_pe.zSum() + ", IMPROVED: " + improvedSolution);
888
889                        if (!improvedSolution) break; else stat_numLSIterationsReducingCost ++; // end grasp iteration 
890                }
891                
892                return current_xrs;
893        }
894        
895        private DoubleMatrix1D computePe_11 (DoubleMatrix1D x_r , DoubleMatrix1D x2_r)
896        {
897                DoubleMatrix1D p_e = A_er.zMult(x_r, null);
898                DoubleMatrix1D backup_occup_e = A_er.zMult(x2_r, null);
899                for (int e = 0 ; e < Efm ; e ++) p_e.set(e , p_e.get(e) + backup_occup_e.get(e));
900                p_e.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil);
901                return p_e;
902        }
903        private DoubleMatrix1D computePe_restoration (DoubleMatrix2D x_rs)
904        {
905                DoubleMatrix2D p_es = A_er.zMult(x_rs, null); DoubleMatrix1D p_e = DoubleFactory1D.dense.make(Efm); 
906                for (int e = 0; e < Efm ; e ++) p_e.set(e , p_es.viewRow(e).getMaxLocation() [0]); // the maximum per state s
907                p_e.assign(DoubleFunctions.div(tcfa_linkCapacity_numCirc.getInt ())).assign(DoubleFunctions.ceil);
908                return p_e;
909        }
910        
911        private double eTime (long init) { return 1E-9*(System.nanoTime() - init); }
912        
913        private static int [] getIndexes (Collection<? extends NetworkElement> col) { final int [] res = new int [col.size ()]; int cont = 0; for (NetworkElement el : col) res [cont++] = el.getIndex (); return res; }
914        
915        
916}