001/*
002 * Copyright 2006 - 2013
003 *     Stefan Balev     <stefan.balev@graphstream-project.org>
004 *     Julien Baudry    <julien.baudry@graphstream-project.org>
005 *     Antoine Dutot    <antoine.dutot@graphstream-project.org>
006 *     Yoann Pigné      <yoann.pigne@graphstream-project.org>
007 *     Guilhelm Savin   <guilhelm.savin@graphstream-project.org>
008 * 
009 * This file is part of GraphStream <http://graphstream-project.org>.
010 * 
011 * GraphStream is a library whose purpose is to handle static or dynamic
012 * graph, create them from scratch, file or any source and display them.
013 * 
014 * This program is free software distributed under the terms of two licenses, the
015 * CeCILL-C license that fits European law, and the GNU Lesser General Public
016 * License. You can  use, modify and/ or redistribute the software under the terms
017 * of the CeCILL-C license as circulated by CEA, CNRS and INRIA at the following
018 * URL <http://www.cecill.info> or under the terms of the GNU LGPL as published by
019 * the Free Software Foundation, either version 3 of the License, or (at your
020 * option) any later version.
021 * 
022 * This program is distributed in the hope that it will be useful, but WITHOUT ANY
023 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
024 * PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
025 * 
026 * You should have received a copy of the GNU Lesser General Public License
027 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
028 * 
029 * The fact that you are presently reading this means that you have had
030 * knowledge of the CeCILL-C and LGPL licenses and that you accept their terms.
031 */
032package org.graphstream.algorithm.networksimplex;
033
034import java.io.PrintStream;
035import java.util.HashMap;
036import java.util.HashSet;
037import java.util.Map;
038import java.util.Set;
039
040import org.graphstream.algorithm.DynamicAlgorithm;
041import org.graphstream.graph.Edge;
042import org.graphstream.graph.Graph;
043import org.graphstream.graph.Node;
044import org.graphstream.stream.SinkAdapter;
045
046// XXX Work in progress
047
048/**
049 * <h3>Minimum cost flow problem</h3>
050 * 
051 * <p>
052 * Network simplex method is an algorithm that solves the minimum cost flow
053 * (MCF) problem for an oriented graph.
054 * </p>
055 * 
056 * <p>
057 * The MCF problem can be stated as follows. Each node has associated number
058 * <i>supply</i> representing the available supply of or demand for flow at that
059 * node. If <i>supply</i> is positive, the node is a supply node, if
060 * <i>supply</i> is negative, the node is a demand node and if <i>supply</i> is
061 * zero, the node is a transshipment node. Each arc has associated
062 * <i>capacity</i> (possibly infinite) and <i>cost</i> per unit flow. The MCF
063 * problem is to send the required flows from supply nodes to demand nodes at
064 * minimum cost, respecting the capacities of the arcs. Note that if the sum of
065 * <i>supply</i> attributes of all nodes is nonzero, the problem is infeasible.
066 * </p>
067 * 
068 * <p>
069 * MCF framework can be used to model a broad variety of network problems,
070 * including matching, shortest path, transportation, etc. For example, if we
071 * want to find the shortest paths from a source to all other nodes in a graph
072 * with <i>n</i> nodes, we can set the <i>supply</i> to <i>n</i>-1 for the
073 * source and to -1 for all other nodes, set <i>capacity</i> to <i>n</i>-1 and
074 * <i>cost</i> to the weight for each arc. The solution of the MCF problem with
075 * these particular settings will be minimum cost unit flow from the source to
076 * all other nodes passing by the shortest paths.
077 * </p>
078 * 
079 * <h3>Problem data</h3>
080 * 
081 * <p>
082 * The user of this class must store the problem data as attributes of the nodes
083 * and the edges of the graph as described below. The names of these attributes
084 * are specified in the constructor
085 * {@link #NetworkSimplex(String, String, String)}. For efficiency reasons all
086 * the data are supposed to be integer. If some of the attributes are real,
087 * their fractional part is ignored. To avoid loss of precision, the user must
088 * scale her data properly if they are real.
089 * </p>
090 * 
091 * <p>
092 * An attribute called {@code supplyName} is used to store the supply (or demand
093 * if negative) of each node. If a node has not an attribute with this name or
094 * if the value of this attribute is not numeric, the node supply is considered
095 * as zero (transshipment node).
096 * </p>
097 * 
098 * <p>
099 * An attribute called {@code capacityName} is used to store the capacity of
100 * each edge. If an edge has not an attribute with this name, or if the value of
101 * this attribute is negative or not numeric, the edge capacity is considered as
102 * infinite.
103 * </p>
104 * 
105 * <p>
106 * An attribute called {@code costName} is used to store the cost per unit flow
107 * of each edge. If an edge has not an attribute with this name or if the value
108 * of this attribute is not numeric, the cost per unit flow of the edge is
109 * considered 1.
110 * </p>
111 * 
112 * <p>
113 * The flow on a directed edge is always from its source node to its target
114 * node. Each undirected edge is considered as a couple of directed edges with
115 * the same capacity and cost per unit flow. In other words, there are possibly
116 * two independent flows on each undirected edge.
117 * </p>
118 * 
119 * <h3>Solutions</h3>
120 * 
121 * TODO
122 * 
123 * <h3>Visualization</h3>
124 * 
125 * TODO
126 * 
127 * @author Stefan Balev
128 */
129public class NetworkSimplex extends SinkAdapter implements DynamicAlgorithm {
130
131        /**
132         * The algorithm maintains some internal data whose names start with this
133         * prefix. The graph <i>must not</i> have any edges whose IDs start with
134         * this prefix.
135         */
136        public static final String PREFIX = "__NS_";
137
138        /**
139         * Used as capacity value for uncapacitated arcs
140         */
141        protected static final int INFINITE_CAPACITY = -1;
142
143        /**
144         * Pricing strategy used at each iteration of the algorithm. Only two simple
145         * strategies are implemented for the moment, more are to come.
146         */
147        public static enum PricingStrategy {
148                /**
149                 * When using this strategy an iteration is faster, but the number of
150                 * iterations is generally bigger
151                 */
152                FIRST_NEGATIVE,
153
154                /**
155                 * When using this strategy an iteration is slower, but the number of
156                 * iterations is generally smaller
157                 */
158                MOST_NEGATIVE
159        }
160
161        /**
162         * The status of the current solution.
163         */
164        public static enum SolutionStatus {
165                /**
166                 * The current solution is outdated. This is the value when the graph
167                 * has changed since the last call of {@link NetworkSimplex#compute()}
168                 */
169                UNDEFINED,
170
171                /**
172                 * The problem is feasible and bounded. The current solution is up to
173                 * date and optimal.
174                 */
175                OPTIMAL,
176
177                /**
178                 * The problem is infeasible, some of the supply/demand constraints
179                 * cannot be satisfied. The current solution is up to date.
180                 */
181                INFEASIBLE,
182
183                /**
184                 * The problem is unbounded (this happens when the graph contains an
185                 * uncapacitated negative cost cycle). The current solution is up to
186                 * date.
187                 */
188                UNBOUNDED
189        }
190
191        /**
192         * Name of the attribute used to store the supply of each node
193         */
194        protected String supplyName;
195
196        /**
197         * Name of the attribute used to store the capacity of each arc
198         */
199        protected String capacityName;
200
201        /**
202         * Name of the attribute used to store the cost per unit flow of each arc
203         */
204        protected String costName;
205
206        /**
207         * Current pricing strategy
208         */
209        protected PricingStrategy pricingStrategy = PricingStrategy.MOST_NEGATIVE;
210
211        /**
212         * A reference to the original graph
213         */
214        protected Graph graph;
215
216        /**
217         * Stores the nodes
218         */
219        protected Map<String, NSNode> nodes;
220
221        /**
222         * Stores the arcs
223         */
224        protected Map<String, NSArc> arcs;
225
226        /**
227         * Stores the non basic arcs.
228         */
229        protected Set<NSArc> nonBasicArcs;
230
231        /**
232         * Artificial root;
233         */
234        protected NSNode root;
235
236        /**
237         * The objective value of the current solution
238         */
239        protected BigMNumber objectiveValue;
240
241        /**
242         * The status of the current BFS
243         */
244        protected SolutionStatus solutionStatus = SolutionStatus.UNDEFINED;
245
246        /**
247         * Entering arc for the next pivot. Set to {@code null} if no candidates.
248         * Set in {@link #selectEnteringArc()}
249         */
250        protected NSArc enteringArc;
251
252        /**
253         * The nearest common predecessor of the extremities of {@link #enteringArc}
254         * . Set in {@link #findJoinNode()}
255         */
256        protected NSNode join;
257
258        /**
259         * The first node of {@link #enteringArc} when traversing it in the
260         * direction of the cycle. Set in {@link #selectLeavingArc()}.
261         */
262        protected NSNode first;
263
264        /**
265         * The second node of {@link #enteringArc} when traversing it in the
266         * direction of the cycle. Set in {@link #selectLeavingArc()}.
267         */
268        protected NSNode second;
269
270        /**
271         * The root of the subtree detached from the BFS tree when removing
272         * {@link #leavingArc}. Set in {@link #selectLeavingArc()}.
273         */
274        protected NSNode oldSubtreeRoot;
275
276        /**
277         * The node of the detached subtree which must be re-attached to the BFS
278         * tree using {@link #enteringArc}. Set in {@link #selectLeavingArc()}.
279         */
280        protected NSNode newSubtreeRoot;
281
282        /**
283         * Maximum allowed flow change on the cycle formed by adding
284         * {@link #enteringArc} to the BFS tree. If infinite, the problem is
285         * unbounded. Set in {@link #selectLeavingArc()}.
286         */
287        protected BigMNumber cycleFlowChange;
288
289        /**
290         * Leaving arc for the next pivot. Set in {@link #selectLeavingArc()}.
291         */
292        protected NSArc leavingArc;
293
294        /**
295         * Working variable, used to avoid frequent instantiations of local
296         * variables.
297         */
298        protected BigMNumber work1;
299
300        /**
301         * Working variable, used to avoid frequent instantiations of local
302         * variables.
303         */
304        protected BigMNumber work2;
305
306        /**
307         * If this delay is positive, sleeps at the end of each pivot and updates UI
308         * classes
309         */
310        protected long animationDelay = 0;
311
312        /**
313         * True if pivot is called from sink method. In this case pivot does not
314         * sleep
315         */
316        protected boolean fromSink = false;
317
318        /**
319         * Log frequency. Logs each logFreq-th iteration
320         */
321        protected int logFreq = 0;
322
323        /**
324         * Log stream
325         */
326        protected PrintStream log = System.err;
327
328        /**
329         * Creates a network simplex instance specifying attribute names to be used.
330         * Use {@link #init(Graph)} to assign a graph to this instance.
331         * 
332         * @param supplyName
333         *            Name of the attribute used to store the supply of each node.
334         * @param capaciyName
335         *            Name of the attribute used to store the capacity of each edge.
336         * @param costName
337         *            Name of the attribute used to store the cost of each edge.
338         */
339        public NetworkSimplex(String supplyName, String capaciyName, String costName) {
340                this.supplyName = supplyName;
341                this.capacityName = capaciyName;
342                this.costName = costName;
343
344                objectiveValue = new BigMNumber();
345                cycleFlowChange = new BigMNumber();
346                work1 = new BigMNumber();
347                work2 = new BigMNumber();
348        }
349
350        // simplex initialization
351
352        /**
353         * Creates copies of all graph arcs and edges. Instantiates and fills
354         * {@link #nodes} and {@link #arcs}.
355         */
356        protected void cloneGraph() {
357                nodes = new HashMap<String, NSNode>(4 * graph.getNodeCount() / 3 + 2);
358                for (Node node : graph) {
359                        NSNode copy = new NSNode(node);
360                        nodes.put(copy.id, copy);
361                }
362
363                int arcCount = graph.getEdgeCount();
364                for (Edge edge : graph.getEachEdge())
365                        if (!edge.isDirected())
366                                arcCount++;
367                arcs = new HashMap<String, NSArc>(4 * arcCount / 3 + 1);
368                for (Edge edge : graph.getEachEdge()) {
369                        NSArc copy = new NSArc(edge, true);
370                        arcs.put(copy.id, copy);
371                        if (!edge.isDirected()) {
372                                copy = new NSArc(edge, false);
373                                arcs.put(copy.id, copy);
374                        }
375                }
376        }
377
378        /**
379         * Creates artificial root and arcs and sets up the initial BFS
380         */
381        protected void createInitialBFS() {
382                nonBasicArcs = new HashSet<NSArc>(4 * arcs.size() / 3 + 1);
383                for (NSArc arc : arcs.values()) {
384                        arc.flow = 0;
385                        arc.status = ArcStatus.NONBASIC_LOWER;
386                        nonBasicArcs.add(arc);
387                        if (animationDelay > 0)
388                                arc.setUIClass();
389                }
390
391                root = new NSNode();
392                root.id = PREFIX + "ROOT";
393                root.potential.set(0);
394                root.parent = root;
395                root.thread = root;
396                root.depth = 0;
397                root.supply = 0;
398                root.artificialArc = null;
399
400                objectiveValue.set(0);
401                for (NSNode node : nodes.values())
402                        node.createArtificialArc();
403                solutionStatus = SolutionStatus.UNDEFINED;
404        }
405
406        // Simplex machinery
407
408        /**
409         * "First negative" pricing strategy
410         */
411        protected void selectEnteringArcFirstNegative() {
412                enteringArc = null;
413                BigMNumber reducedCost = work1;
414
415                for (NSArc arc : nonBasicArcs) {
416                        arc.computeReducedCost(reducedCost);
417                        if (reducedCost.isNegative()) {
418                                enteringArc = arc;
419                                return;
420                        }
421                }
422                
423                // Skip the artificial arcs if the objective value is finite
424                if (!objectiveValue.isInfinite())
425                        return;
426
427                for (NSNode node : nodes.values()) {
428                        NSArc arc = node.artificialArc;
429                        if (arc.status == ArcStatus.NONBASIC_LOWER) {
430                                arc.computeReducedCost(reducedCost);
431                                if (reducedCost.isNegative()) {
432                                        enteringArc = arc;
433                                        return;
434                                }
435                        }
436                }
437        }
438
439        /**
440         * "Most negative" pricing strategy
441         */
442        protected void selectEnteringArcMostNegative() {
443                enteringArc = null;
444                BigMNumber reducedCost = work1;
445                BigMNumber bestReducedCost = work2;
446                bestReducedCost.set(0);
447
448                for (NSArc arc : nonBasicArcs) {
449                        arc.computeReducedCost(reducedCost);
450                        if (reducedCost.compareTo(bestReducedCost) < 0) {
451                                bestReducedCost.set(reducedCost);
452                                enteringArc = arc;
453                        }
454                }
455                if (enteringArc != null)
456                        return;
457                
458                // Skip the artificial arcs if the objective value is finite
459                if (!objectiveValue.isInfinite())
460                        return;
461
462                for (NSNode node : nodes.values()) {
463                        NSArc arc = node.artificialArc;
464                        if (arc.status == ArcStatus.NONBASIC_LOWER) {
465                                arc.computeReducedCost(reducedCost);
466                                if (reducedCost.compareTo(bestReducedCost) < 0) {
467                                        bestReducedCost = reducedCost;
468                                        enteringArc = arc;
469                                }
470                        }
471                }
472        }
473
474        /**
475         * Selects entering arc among the candidates (non-basic arcs with negative
476         * reduced costs). Puts the selected candidate in {@link #enteringArc}. If
477         * there are no candidates, {@link #enteringArc} is set to {@code null}.
478         */
479        protected void selectEnteringArc() {
480                switch (pricingStrategy) {
481                case FIRST_NEGATIVE:
482                        selectEnteringArcFirstNegative();
483                        break;
484                case MOST_NEGATIVE:
485                        selectEnteringArcMostNegative();
486                        break;
487                }
488        }
489
490        /**
491         * Finds the nearest common predecessor of the two nodes of
492         * {@link #enteringArc}. Puts it in {@link #join}.
493         */
494        protected void findJoinNode() {
495                NSNode i = enteringArc.source;
496                NSNode j = enteringArc.target;
497                while (i.depth > j.depth)
498                        i = i.parent;
499                while (j.depth > i.depth)
500                        j = j.parent;
501                while (i != j) {
502                        i = i.parent;
503                        j = j.parent;
504                }
505                join = i;
506        }
507
508        /**
509         * Selects the leaving arc, that is the arc belonging to the cycle with
510         * minimum allowed flow change. Maintains strongly feasible basis: if there
511         * are more than one candidates, selects the last visited when traversing
512         * the cycle in its direction starting from {@link join}. Sets up
513         * {@link #first}, {@link second}, {@link #cycleFlowChange},
514         * {@link #oldSubtreeRoot}, {@link #newSubtreeRoot} and {@link #leavingArc}.
515         */
516        protected void selectLeavingArc() {
517                findJoinNode();
518                if (enteringArc.status == ArcStatus.NONBASIC_LOWER) {
519                        first = enteringArc.source;
520                        second = enteringArc.target;
521                } else {
522                        first = enteringArc.target;
523                        second = enteringArc.source;
524                }
525
526                enteringArc.computeAllowedFlowChange(first, cycleFlowChange);
527                leavingArc = enteringArc;
528
529                NSArc arc;
530                BigMNumber arcFlowChange = work1;
531
532                for (NSNode node = second; node != join; node = node.parent) {
533                        arc = node.arcToParent;
534                        arc.computeAllowedFlowChange(node, arcFlowChange);
535                        if (arcFlowChange.compareTo(cycleFlowChange) <= 0) {
536                                cycleFlowChange.set(arcFlowChange);
537                                oldSubtreeRoot = node;
538                                newSubtreeRoot = second;
539                                leavingArc = arc;
540                        }
541                }
542
543                for (NSNode node = first; node != join; node = node.parent) {
544                        arc = node.arcToParent;
545                        arc.computeAllowedFlowChange(node.parent, arcFlowChange);
546                        if (arcFlowChange.compareTo(cycleFlowChange) < 0) {
547                                cycleFlowChange.set(arcFlowChange);
548                                oldSubtreeRoot = node;
549                                newSubtreeRoot = first;
550                                leavingArc = arc;
551                        }
552                }
553        }
554
555        /**
556         * Changes the flows on the arcs belonging to the cycle and updates the
557         * objective value.
558         */
559        protected void changeFlows() {
560                int delta = (int) cycleFlowChange.getSmall();
561                if (delta == 0)
562                        return;
563
564                enteringArc.computeReducedCost(work1);
565                objectiveValue.plusTimes(delta, work1);
566
567                enteringArc.changeFlow(delta, first);
568                for (NSNode node = second; node != join; node = node.parent)
569                        node.arcToParent.changeFlow(delta, node);
570                for (NSNode node = first; node != join; node = node.parent)
571                        node.arcToParent.changeFlow(delta, node.parent);
572        }
573
574        /**
575         * Turns upside-down the part of the tree between the entering arc and the
576         * leaving arc.
577         */
578        protected void updateBFS() {
579                NSNode stopNode = oldSubtreeRoot.parent;
580
581                NSNode currentNode = newSubtreeRoot;
582                NSNode oldParent = currentNode.parent;
583                NSNode newParent = enteringArc.getOpposite(currentNode);
584                NSArc oldArc = currentNode.arcToParent;
585                NSArc newArc = enteringArc;
586                while (currentNode != stopNode) {
587                        currentNode.changeParent(newParent, newArc);
588                        newParent = currentNode;
589                        currentNode = oldParent;
590                        oldParent = currentNode.parent;
591                        newArc = oldArc;
592                        oldArc = currentNode.arcToParent;
593                }
594        }
595
596        /**
597         * Performs a pivot when the entering arc and the leaving arc are known.
598         * Changes the flows on the cycle and updates the BFS.
599         */
600        protected void pivot() {
601                if (animationDelay > 0 && !fromSink)
602                        try {
603                                Thread.sleep(animationDelay);
604                        } catch (InterruptedException e) {
605                        }
606
607                changeFlows();
608                if (enteringArc == leavingArc) {
609                        if (enteringArc.status == ArcStatus.NONBASIC_LOWER)
610                                enteringArc.status = ArcStatus.NONBASIC_UPPER;
611                        else
612                                enteringArc.status = ArcStatus.NONBASIC_LOWER;
613                } else {
614                        enteringArc.status = ArcStatus.BASIC;
615                        nonBasicArcs.remove(enteringArc);
616                        if ((newSubtreeRoot == first && oldSubtreeRoot == leavingArc.target)
617                                        || (newSubtreeRoot == second && oldSubtreeRoot == leavingArc.source))
618                                // The leaving arc is in the direction of the cycle
619                                leavingArc.status = ArcStatus.NONBASIC_UPPER;
620                        else
621                                leavingArc.status = ArcStatus.NONBASIC_LOWER;
622                        if (!leavingArc.isArtificial())
623                                nonBasicArcs.add(leavingArc);
624                        updateBFS();
625                }
626                if (animationDelay > 0) {
627                        enteringArc.setUIClass();
628                        leavingArc.setUIClass();
629                }
630        }
631
632        /**
633         * The main simplex method loop. Selects leaving and entering arc and
634         * performs a pivot. Loops until there are no more candidates or until
635         * absorbing cycle is found.
636         */
637        protected void simplex() {
638                int pivots = 0;
639                if (logFreq > 0) {
640                        log.println("Starting simplex...");
641                        log.printf("%10s%30s%30s%10s%10s%10s%n", "pivot", "entering",
642                                        "leaving", "delta", "cost", "infeas.");
643                }
644                while (true) {
645                        selectEnteringArc();
646                        if (enteringArc == null) {
647                                if (objectiveValue.isInfinite())
648                                        solutionStatus = SolutionStatus.INFEASIBLE;
649                                else
650                                        solutionStatus = SolutionStatus.OPTIMAL;
651                                break;
652                        }
653                        selectLeavingArc();
654                        if (cycleFlowChange.isInfinite()) {
655                                solutionStatus = SolutionStatus.UNBOUNDED;
656                                break;
657                        }
658                        pivot();
659                        pivots++;
660                        if (logFreq > 0 && pivots % logFreq == 0)
661                                log.printf("%10d%30s%30s%10d%10d%10d%n", pivots,
662                                                enteringArc.id, leavingArc.id, cycleFlowChange.small,
663                                                objectiveValue.small, objectiveValue.big);
664                }
665                if (logFreq > 0)
666                        log.printf(
667                                        "Simplex finished (%d pivots). Cost: %d. Status: %s%n%n",
668                                        pivots, objectiveValue.small, solutionStatus);
669        }
670
671        // access and modification of algorithm parameters
672
673        /**
674         * Returns the name of the attribute used to store the supply of each node.
675         * This name is given as constructor parameter and cannot be modified.
676         * 
677         * @return The name of the supply attribute.
678         */
679        public String getSupplyName() {
680                return supplyName;
681        }
682
683        /**
684         * Returns the name of the attribute used to store the capacity of each
685         * edge. This name is given as constructor parameter and cannot be modified.
686         * 
687         * @return The name of the capacity attribute.
688         */
689        public String getCapacityName() {
690                return capacityName;
691        }
692
693        /**
694         * Returns the name of the attribute used to store the cost per unit flow of
695         * each edge. This name is given as constructor parameter and cannot be
696         * modified.
697         * 
698         * @return The name of the cost attribute.
699         */
700        public String getCostName() {
701                return costName;
702        }
703
704        /**
705         * Returns the currently used pricing strategy.
706         * 
707         * @return The pricing strategy
708         */
709        public PricingStrategy getPricingStrategy() {
710                return pricingStrategy;
711        }
712
713        /**
714         * Sets the pricing strategy
715         * 
716         * @param pricingStrategy
717         *            The new pricing strategy
718         */
719        public void setPricingStrategy(PricingStrategy pricingStrategy) {
720                this.pricingStrategy = pricingStrategy;
721        }
722
723        /**
724         * When the animation delay is positive, the algorithm continuously updates
725         * {@code "ui.class"} and {@code "label"} attributes of the edges and the
726         * nodes of the graph and sleeps at the beginning of each simplex pivot.
727         * This feature can be useful for visualizing the algorithm execution. The
728         * user must provide a stylesheet defining the classes of the graph elements
729         * as described in {@link #setUIClasses()}. This feature is disabled by
730         * default.
731         * 
732         * @param millis
733         *            The time in milliseconds to sleep between two simplex pivots.
734         * @see #setUIClasses()
735         */
736        public void setAnimationDelay(long millis) {
737                animationDelay = millis;
738        }
739
740        /**
741         * Returns the graph on which the algorithm is applied. This is the graph
742         * passed in parameter in {@link #init(Graph)}.
743         * 
744         * @return The graph on which the algorithm is applied.
745         */
746        public Graph getGraph() {
747                return graph;
748        }
749
750        /**
751         * Sets the log frequency.
752         * 
753         * If the parameter is positive, outputs information about the algorithm
754         * execution to the log stream.
755         * 
756         * @param pivots
757         *            The log frequency in number of pivots
758         * @see #setLogStream(PrintStream)
759         */
760        public void setLogFrequency(int pivots) {
761                logFreq = pivots;
762        }
763
764        /**
765         * Sets the log stream.
766         * 
767         * Note that the algorithm outputs information about its execution only if
768         * the log frequency is positive. By default the log stream is
769         * {@link System#err}.
770         * 
771         * @param log
772         *            The log stream
773         * @see #setLogFrequency(int)
774         */
775        public void setLogStream(PrintStream log) {
776                this.log = log;
777        }
778
779        /**
780         * Returns the sum of the supplies of all the nodes in the network.
781         * 
782         * The MCF problem has solution only if the problem is balanced, i.e. if the
783         * total supply is equal to the total demand. This method returns the
784         * missing supply (if negative) or demand (if positive) in order to make the
785         * problem balanced. If the returned value is zero, the problem is balanced.
786         * 
787         * @return The network balance
788         */
789        public int getNetworkBalance() {
790                return -root.supply;
791        }
792
793        // solution access methods
794
795        /**
796         * If the current solution is up to date, returns the status of the problem.
797         * Otherwise returns {@link SolutionStatus#UNDEFINED}.
798         * 
799         * 
800         * @return The status of the current solution.
801         * @see SolutionStatus
802         */
803        public SolutionStatus getSolutionStatus() {
804                return solutionStatus;
805        }
806
807        /**
808         * Returns the total cost of the current network flow
809         * 
810         * @return The cost of the flow defined by the current solution
811         */
812        public long getSolutionCost() {
813                return objectiveValue.getSmall();
814        }
815
816        /**
817         * Returns the infeasibility of the current solution.
818         * 
819         * This is the sum of the absolute values of the infeasibilities of all the
820         * nodes. If the returned value is zero, the current solution is feasible,
821         * i.e. it satisfies the supply constraints of all the nodes.
822         * 
823         * @return The infeasibility of the current solution.
824         * @see #getInfeasibility(Node)
825         */
826        public long getSolutionInfeasibility() {
827                return objectiveValue.big;
828        }
829
830        /**
831         * Returns the infeasibility of a node.
832         * 
833         * Returns the amount of missing outflow (if positive) or inflow (if
834         * negative) of a given node. If the value is zero, the current solution
835         * satisfies the node demand / supply.
836         * 
837         * @param node
838         *            A node
839         * @return The infeasibility of the node
840         */
841        public int getInfeasibility(Node node) {
842                NSArc artificial = nodes.get(node.getId()).artificialArc;
843                return artificial.target == root ? artificial.flow : -artificial.flow;
844        }
845
846        /**
847         * Returns the edge to the parent of a node in the current BFS tree.
848         * 
849         * If the parent of the node is the artificial root, this method returns
850         * {@code null}. When the returned edge is undirected, use
851         * {@link #getStatus(Edge, boolean)} to know which of the both arcs is
852         * basic.
853         * 
854         * @param node
855         *            A node
856         * @return The edge to the parent of the node in the BFS tree
857         */
858        public <T extends Edge> T getEdgeFromParent(Node node) {
859                NSArc arc = nodes.get(node.getId()).arcToParent;
860                if (arc.isArtificial())
861                        return null;
862                return graph.getEdge(arc.getOriginalId());
863        }
864
865        /**
866         * Returns the parent of a node in the current BFS tree.
867         * 
868         * If the parent of the node is the artificial root, returns {@code null}.
869         * 
870         * @param node
871         *            A node
872         * @return The parent of a node in the BFS tree
873         */
874        public <T extends Node> T getParent(Node node) {
875                NSNode nsNode = nodes.get(node.getId());
876                if (nsNode == root)
877                        return null;
878                return graph.getNode(nsNode.parent.id);
879        }
880
881        /**
882         * Returns the flow on an edge.
883         * 
884         * If {@code sameDirection} is true, returns the flow from the source to the
885         * target of the edge, otherwise returns the flow from the target to the
886         * source of the edge. Note that for directed edges the flow can only pass
887         * from the source node to the target node. For undirected edges there may
888         * be independent flows in both directions.
889         * 
890         * @param edge
891         *            An edge
892         * @param sameDirection
893         *            If true, returns the flow from the source to the target.
894         * @return The flow on the edge.
895         */
896        public int getFlow(Edge edge, boolean sameDirection) {
897                if (edge.isDirected())
898                        return sameDirection ? arcs.get(edge.getId()).flow : 0;
899                else
900                        return arcs.get((sameDirection ? "" : PREFIX + "REVERSE_")
901                                        + edge.getId()).flow;
902        }
903
904        /**
905         * Returns the flow on an edge from its source node to its target node.
906         * 
907         * The same as {@code getFlow(Edge, true)}.
908         * 
909         * @param edge
910         *            An edge
911         * @return The flow on the edge
912         * @see #getFlow(Edge, boolean)
913         */
914        public int getFlow(Edge edge) {
915                return getFlow(edge, true);
916        }
917
918        /**
919         * Returns the status of an edge in the current solution.
920         * 
921         * An edge can be basic, non-basic at zero or non-basic at upper bound. Note
922         * that undirected edges are interpreted as two directed arcs. If
923         * {@code sameDirection} is true, the method returns the status of the arc
924         * from the source to the target of the edge, otherwise it returns the
925         * status of the arc from the target to the source. If the edge is directed
926         * and {@code sameDirection} is false, returns {@code null}.
927         * 
928         * @param edge
929         *            An edge
930         * @param sameDirection
931         *            If true, returns the status of the arc from the source to the
932         *            target.
933         * @return The status of the edge
934         */
935        public ArcStatus getStatus(Edge edge, boolean sameDirection) {
936                if (edge.isDirected())
937                        return sameDirection ? arcs.get(edge.getId()).status : null;
938                else
939                        return arcs.get((sameDirection ? "" : PREFIX + "REVERSE_")
940                                        + edge.getId()).status;
941        }
942
943        /**
944         * Returns the status of an edge in the current solution.
945         * 
946         * The same as {@code getStatus(edge, true)}.
947         * 
948         * @param edge
949         *            An edge
950         * @return The status of the edge
951         * @see #getStatus(Edge, boolean)
952         */
953        public ArcStatus getStatus(Edge edge) {
954                return getStatus(edge, true);
955        }
956
957        /**
958         * This method can be used to visualize the current solution.
959         * 
960         * <p>
961         * It sets the attributes {@code "label"} and {@code "ui.class"} of the
962         * nodes and the edges of the graph depending on the current solution. The
963         * labels of the nodes are set to their balance (see
964         * {@link #getInfeasibility(Node)}). The labels of the edges are set to the
965         * flow passing through them. The {@code "ui.class"} attribute of the nodes
966         * is set to one of {@code "supply_balanced"}, {@code "supply_unbalanced"},
967         * {@code "demand_balanced"}, {@code "demand_unbalanced"},
968         * {@code "trans_balanced"} or {@code "trans_unbalanced"} depending on the
969         * node type and the node status. The {@code "ui.class"} attribute of the
970         * edges is set to one of {@code "basic"}, {@code "nonbasic_lower"} or
971         * {@code "nonbasic_upper"} according to their status (see
972         * {@link #getStatus(Edge)}.
973         * </p>
974         * <p>
975         * The user must provide a stylesheet defining the visual appearance for
976         * each of these node and edge classes. Note that if the animation delay is
977         * positive (see {@link #setAnimationDelay(long)}), there is no need to call
978         * this method, because in this case the labels and the UI classes of the
979         * graph elements are set and updated during the algorithm execution.
980         * </p>
981         * <p>
982         * Note that in the case of undirected edges the label and the UI class are
983         * set according to the status of one of the corresponding arcs (not
984         * specified which one).
985         * </p>
986         */
987        public void setUIClasses() {
988                for (Node node : graph)
989                        nodes.get(node.getId()).artificialArc.setUIClass();
990                for (Edge edge : graph.getEachEdge()) {
991                        NSArc arc = arcs.get(edge.getId());
992                        if (!edge.isDirected() && arc.status != ArcStatus.BASIC)
993                                arc = arcs.get(PREFIX + "REVERSE_" + edge.getId());
994                        arc.setUIClass();
995                }
996        }
997
998        // DynamicAlgorithm methods
999
1000        public void init(Graph graph) {
1001                this.graph = graph;
1002                cloneGraph();
1003                createInitialBFS();
1004                graph.addSink(this);
1005        }
1006
1007        public void compute() {
1008                fromSink = false;
1009                if (solutionStatus == SolutionStatus.UNDEFINED)
1010                        simplex();
1011        }
1012
1013        public void terminate() {
1014                graph.removeSink(this);
1015                solutionStatus = SolutionStatus.UNDEFINED;
1016        }
1017
1018        // Sink methods
1019
1020        @Override
1021        public void edgeAttributeAdded(String sourceId, long timeId, String edgeId,
1022                        String attribute, Object value) {
1023                if (attribute.equals(costName)) {
1024                        double v = objectToDouble(value);
1025                        if (Double.isNaN(v))
1026                                v = 1;
1027
1028                        NSArc arc = arcs.get(edgeId);
1029                        work1.set((int) v);
1030                        changeCost(arc, work1);
1031
1032                        arc = arcs.get(PREFIX + "REVERSE_" + edgeId);
1033                        if (arc != null) {
1034                                work1.set((int) v);
1035                                changeCost(arc, work1);
1036                        }
1037                } else if (attribute.equals(capacityName)) {
1038                        double v = objectToDouble(value);
1039                        if (Double.isNaN(v) || v < 0)
1040                                v = INFINITE_CAPACITY;
1041                        NSArc arc = arcs.get(edgeId);
1042                        changeCapacity(arc, (int) v);
1043                        arc = arcs.get(PREFIX + "REVERSE_" + edgeId);
1044                        if (arc != null)
1045                                changeCapacity(arc, (int) v);
1046                }
1047        }
1048
1049        @Override
1050        public void edgeAttributeChanged(String sourceId, long timeId,
1051                        String edgeId, String attribute, Object oldValue, Object newValue) {
1052                edgeAttributeAdded(sourceId, timeId, edgeId, attribute, newValue);
1053        }
1054
1055        @Override
1056        public void edgeAttributeRemoved(String sourceId, long timeId,
1057                        String edgeId, String attribute) {
1058                edgeAttributeAdded(sourceId, timeId, edgeId, attribute, null);
1059        }
1060
1061        @Override
1062        public void nodeAttributeAdded(String sourceId, long timeId, String nodeId,
1063                        String attribute, Object value) {
1064                if (attribute.equals(supplyName)) {
1065                        double v = objectToDouble(value);
1066                        if (Double.isNaN(v))
1067                                v = 0;
1068
1069                        NSNode node = nodes.get(nodeId);
1070                        changeSupply(node, (int) v);
1071                }
1072        }
1073
1074        @Override
1075        public void nodeAttributeChanged(String sourceId, long timeId,
1076                        String nodeId, String attribute, Object oldValue, Object newValue) {
1077                nodeAttributeAdded(sourceId, timeId, nodeId, attribute, newValue);
1078        }
1079
1080        @Override
1081        public void nodeAttributeRemoved(String sourceId, long timeId,
1082                        String nodeId, String attribute) {
1083                nodeAttributeAdded(sourceId, timeId, nodeId, attribute, null);
1084        }
1085
1086        @Override
1087        public void edgeAdded(String sourceId, long timeId, String edgeId,
1088                        String fromNodeId, String toNodeId, boolean directed) {
1089                NSArc arc = new NSArc(graph.getEdge(edgeId), true);
1090                addArc(arc);
1091                if (!directed) {
1092                        arc = new NSArc(graph.getEdge(edgeId), false);
1093                        addArc(arc);
1094                }
1095        }
1096
1097        @Override
1098        public void edgeRemoved(String sourceId, long timeId, String edgeId) {
1099                NSArc arc = arcs.get(edgeId);
1100                removeArc(arc);
1101                arc = arcs.get(PREFIX + "REVERSE_" + edgeId);
1102                if (arc != null)
1103                        removeArc(arc);
1104        }
1105
1106        @Override
1107        public void nodeAdded(String sourceId, long timeId, String nodeId) {
1108                addNode(new NSNode(graph.getNode(nodeId)));
1109        }
1110
1111        @Override
1112        public void nodeRemoved(String sourceId, long timeId, String nodeId) {
1113                removeNode(nodes.get(nodeId));
1114        }
1115
1116        @Override
1117        public void graphCleared(String sourceId, long timeId) {
1118                clearGraph();
1119        }
1120
1121        // helpers for the sink
1122
1123        /**
1124         * Utility method trying to convert object to double in the same way as
1125         * {@link org.graphstream.graph.implementations.AbstractElement#getNumber(String)}
1126         * does.
1127         * 
1128         * @param o
1129         *            The object to be converted
1130         * @return The numeric value of the object or NaN
1131         */
1132        protected static double objectToDouble(Object o) {
1133                if (o != null) {
1134                        if (o instanceof Number)
1135                                return ((Number) o).doubleValue();
1136
1137                        if (o instanceof String) {
1138                                try {
1139                                        return Double.parseDouble((String) o);
1140                                } catch (NumberFormatException e) {
1141                                }
1142                        }
1143                }
1144                return Double.NaN;
1145        }
1146
1147        /**
1148         * Changes the cost of an arc
1149         * 
1150         * @param arc
1151         *            The arc that changes cost
1152         * @param newCost
1153         *            The new cost
1154         */
1155        protected void changeCost(NSArc arc, BigMNumber newCost) {
1156                if (arc.cost.compareTo(newCost) == 0)
1157                        return;
1158                objectiveValue.plusTimes(-arc.flow, arc.cost);
1159                arc.cost.set(newCost);
1160                objectiveValue.plusTimes(arc.flow, arc.cost);
1161
1162                if (arc.status == ArcStatus.BASIC) {
1163                        NSNode subtreeRoot = arc.source.arcToParent == arc ? arc.source
1164                                        : arc.target;
1165                        subtreeRoot.computePotential();
1166                        for (NSNode node = subtreeRoot.thread; node.depth > subtreeRoot.depth; node = node.thread)
1167                                node.computePotential();
1168                        solutionStatus = SolutionStatus.UNDEFINED;
1169                } else {
1170                        arc.computeReducedCost(work1);
1171                        if (work1.isNegative())
1172                                solutionStatus = SolutionStatus.UNDEFINED;
1173                }
1174        }
1175
1176        protected void changeSupply(NSNode node, int newSupply) {
1177                if (node.supply == newSupply)
1178                        return;
1179                NSArc artificial = node.artificialArc;
1180                // enter the artificial arc in the tree if not there
1181                if (artificial.status == ArcStatus.NONBASIC_LOWER) {
1182                        enteringArc = artificial;
1183                        selectLeavingArc();
1184                        // if we are in infinite cycle, switch the direction
1185                        if (cycleFlowChange.isInfinite()) {
1186                                artificial.switchDirection();
1187                                selectLeavingArc();
1188                        }
1189                        fromSink = true;
1190                        pivot();
1191                }
1192                // now the artificial arc is basic and we can change its flow
1193                objectiveValue.plusTimes(-artificial.flow, artificial.cost);
1194                int delta = newSupply - node.supply;
1195                node.supply = newSupply;
1196                root.supply -= delta;
1197                if (node == artificial.source) {
1198                        artificial.flow += delta;
1199                } else {
1200                        artificial.flow -= delta;
1201                }
1202                if (artificial.flow < 0)
1203                        artificial.switchDirection();
1204
1205                objectiveValue.plusTimes(artificial.flow, artificial.cost);
1206                solutionStatus = SolutionStatus.UNDEFINED;
1207
1208                if (animationDelay > 0)
1209                        artificial.setUIClass();
1210        }
1211
1212        protected void changeCapacity(NSArc arc, int newCapacity) {
1213                if (arc.capacity == newCapacity)
1214                        return;
1215                if (arc.status == ArcStatus.NONBASIC_LOWER) {
1216                        arc.capacity = newCapacity;
1217                        return;
1218                }
1219                if (arc.status == ArcStatus.NONBASIC_UPPER) {
1220                        enteringArc = arc;
1221                        selectLeavingArc();
1222                        fromSink = true;
1223                        pivot();
1224                        solutionStatus = SolutionStatus.UNDEFINED;
1225                }
1226                // now the arc is basic ...
1227                if (newCapacity == INFINITE_CAPACITY || arc.flow <= newCapacity) {
1228                        arc.capacity = newCapacity;
1229                        return;
1230                }
1231                // ... and the flow on it is greater than its new capacity
1232                int delta = arc.flow - newCapacity;
1233                arc.flow = arc.capacity = newCapacity;
1234                objectiveValue.plusTimes(-delta, arc.cost);
1235                arc.source.supply -= delta;
1236                arc.target.supply += delta;
1237                if (animationDelay > 0)
1238                        arc.setUIClass();
1239                changeSupply(arc.source, arc.source.supply + delta);
1240                changeSupply(arc.target, arc.target.supply - delta);
1241        }
1242
1243        protected void addArc(NSArc arc) {
1244                arc.flow = 0;
1245                arc.status = ArcStatus.NONBASIC_LOWER;
1246                arcs.put(arc.id, arc);
1247                nonBasicArcs.add(arc);
1248                arc.computeReducedCost(work1);
1249                if (work1.isNegative())
1250                        solutionStatus = SolutionStatus.UNDEFINED;
1251                if (animationDelay > 0)
1252                        arc.setUIClass();
1253        }
1254
1255        protected void removeArc(NSArc arc) {
1256                changeCapacity(arc, 0);
1257                if (arc.status == ArcStatus.BASIC) {
1258                        NSNode node = arc.source.arcToParent == arc ? arc.source
1259                                        : arc.target;
1260                        enteringArc = node.artificialArc;
1261                        if (enteringArc.source == root)
1262                                enteringArc.switchDirection();
1263                        selectLeavingArc();
1264                        fromSink = true;
1265                        pivot();
1266                        solutionStatus = SolutionStatus.UNDEFINED;
1267                }
1268                arcs.remove(arc.id);
1269                nonBasicArcs.remove(arc);
1270        }
1271
1272        protected void addNode(NSNode node) {
1273                nodes.put(node.id, node);
1274                node.createArtificialArc();
1275                solutionStatus = SolutionStatus.UNDEFINED;
1276        }
1277
1278        protected void removeNode(NSNode node) {
1279                node.previousInThread().thread = node.thread;
1280                NSArc artificial = node.arcToParent;
1281                objectiveValue.plusTimes(-artificial.flow, artificial.cost);
1282                root.supply += node.supply;
1283                nodes.remove(node.id);
1284                solutionStatus = SolutionStatus.UNDEFINED;
1285        }
1286
1287        protected void clearGraph() {
1288                nodes.clear();
1289                arcs.clear();
1290                nonBasicArcs.clear();
1291                root.thread = root;
1292                root.supply = 0;
1293                objectiveValue.set(0);
1294                solutionStatus = SolutionStatus.OPTIMAL;
1295        }
1296
1297        /**
1298         * Internal representation of the graph nodes. Stores node ids, supplies and
1299         * potentials. Maintains BFS tree using PTD (parent, thread, depth) data
1300         * structure.
1301         */
1302        protected class NSNode {
1303
1304                /**
1305                 * Node id. The same as in the original graph. Special id for the
1306                 * artificial root.
1307                 */
1308                String id;
1309
1310                /**
1311                 * Node supply (or demand if negative). Problem data retrieved from the
1312                 * original graph.
1313                 */
1314                int supply;
1315
1316                /**
1317                 * Node potential
1318                 */
1319                BigMNumber potential;
1320
1321                /**
1322                 * Parent in the BFS tree
1323                 */
1324                NSNode parent;
1325
1326                /**
1327                 * The next node in the preorder traversal of the BFS tree
1328                 */
1329                NSNode thread;
1330
1331                /**
1332                 * The depth in the BFS tree
1333                 */
1334                int depth;
1335
1336                /**
1337                 * The arc connecting this node to its parent in the BFS tree
1338                 */
1339                NSArc arcToParent;
1340
1341                /**
1342                 * The artificial arc associated to this node
1343                 */
1344                NSArc artificialArc;
1345
1346                /**
1347                 * Creates a copy of a node.
1348                 * 
1349                 * @param node
1350                 *            a node of the original graph
1351                 */
1352                NSNode(Node node) {
1353                        id = node.getId();
1354                        double v = node.getNumber(supplyName);
1355                        if (Double.isNaN(v))
1356                                v = 0;
1357                        supply = (int) v;
1358                        potential = new BigMNumber();
1359                }
1360
1361                /**
1362                 * Default constructor.
1363                 */
1364                NSNode() {
1365                        potential = new BigMNumber();
1366                }
1367
1368                /**
1369                 * Creates the artificial arc corresponding to this node and puts it in
1370                 * the BFS
1371                 */
1372                void createArtificialArc() {
1373                        artificialArc = new NSArc();
1374                        artificialArc.id = id;
1375                        artificialArc.capacity = INFINITE_CAPACITY;
1376                        artificialArc.cost.set(0, 1);
1377                        artificialArc.status = ArcStatus.BASIC;
1378                        if (supply > 0) {
1379                                artificialArc.source = this;
1380                                artificialArc.target = root;
1381                                artificialArc.flow = supply;
1382                        } else {
1383                                artificialArc.source = root;
1384                                artificialArc.target = this;
1385                                artificialArc.flow = -supply;
1386                        }
1387
1388                        parent = root;
1389                        thread = root.thread;
1390                        root.thread = this;
1391                        depth = 1;
1392                        arcToParent = artificialArc;
1393                        computePotential();
1394
1395                        root.supply -= supply;
1396                        objectiveValue.plusTimes(artificialArc.flow, artificialArc.cost);
1397
1398                        if (animationDelay > 0)
1399                                artificialArc.setUIClass();
1400                }
1401
1402                /**
1403                 * Finds the previous node in the preorder traversal of the BFS tree
1404                 * 
1405                 * @return the previous node in the thread
1406                 */
1407                NSNode previousInThread() {
1408                        NSNode node;
1409                        for (node = parent; node.thread != this; node = node.thread)
1410                                ;
1411                        return node;
1412                }
1413
1414                /**
1415                 * Finds the rightmost node of the subtree of this node when following
1416                 * the thread.
1417                 * 
1418                 * @return The last successor of this node
1419                 */
1420                NSNode lastSuccessor() {
1421                        NSNode node;
1422                        for (node = this; node.thread.depth > depth; node = node.thread)
1423                                ;
1424                        return node;
1425                }
1426
1427                /**
1428                 * Computes the potential of this node knowing the potential of its
1429                 * father
1430                 */
1431                void computePotential() {
1432                        potential.set(parent.potential);
1433                        if (arcToParent.source == this)
1434                                potential.plus(arcToParent.cost);
1435                        else
1436                                potential.minus(arcToParent.cost);
1437                }
1438
1439                /**
1440                 * Changes the parent of this node. Updates PTD structure and node
1441                 * potentials.
1442                 * 
1443                 * @param newParent
1444                 *            the new parent
1445                 * @param newArcToParent
1446                 *            the arc to the new parent
1447                 */
1448                void changeParent(NSNode newParent, NSArc newArcToParent) {
1449                        NSNode pred = previousInThread();
1450                        NSNode succ = lastSuccessor();
1451
1452                        pred.thread = succ.thread;
1453                        succ.thread = newParent.thread;
1454                        newParent.thread = this;
1455
1456                        parent = newParent;
1457                        arcToParent = newArcToParent;
1458
1459                        for (NSNode node = this; node != succ.thread; node = node.thread) {
1460                                node.depth = node.parent.depth + 1;
1461                                node.computePotential();
1462                        }
1463                }
1464        }
1465
1466        /**
1467         * Arc status
1468         */
1469        public static enum ArcStatus {
1470                /**
1471                 * Basic arc
1472                 */
1473                BASIC,
1474                /**
1475                 * Non basic arc with zero flow
1476                 */
1477                NONBASIC_LOWER,
1478                /**
1479                 * Non basic saturated arc
1480                 */
1481                NONBASIC_UPPER;
1482        }
1483
1484        /**
1485         * Internal representation of the graph arcs. Stores the arc ids,
1486         * capacities, costs, source and target nodes. Maintains BFS information:
1487         * flow and status.
1488         */
1489        protected class NSArc {
1490                /**
1491                 * Arc id. The same as in the original graph. Special ids for the
1492                 * artificial arcs and the arcs doubling undirected edges.
1493                 */
1494                String id;
1495
1496                /**
1497                 * Arc capacity. Problem data retrieved from the original graph.
1498                 * Considered infinite if equal to -1.
1499                 */
1500                int capacity;
1501
1502                /**
1503                 * Arc cost. Problem data retrieved from the original graph. Big M for
1504                 * the artificial arcs.
1505                 */
1506                BigMNumber cost;
1507
1508                /**
1509                 * Source node
1510                 */
1511                NSNode source;
1512
1513                /**
1514                 * Target node
1515                 */
1516                NSNode target;
1517
1518                /**
1519                 * Flow on this arc in the current BFS
1520                 */
1521                int flow;
1522
1523                /**
1524                 * Status of this arc in the current BFS
1525                 */
1526                ArcStatus status;
1527
1528                /**
1529                 * Creates a copy of an edge of the original graph
1530                 * 
1531                 * @param edge
1532                 *            an edge of the original graph
1533                 * @param sameDirection
1534                 *            true if the arc must have the same direction as the
1535                 *            original edge
1536                 */
1537                NSArc(Edge edge, boolean sameDirection) {
1538                        if (edge.getId().startsWith(PREFIX))
1539                                throw new IllegalArgumentException(
1540                                                "Edge ids must not start with " + PREFIX);
1541                        id = (sameDirection ? "" : PREFIX + "REVERSE_") + edge.getId();
1542
1543                        double v = edge.getNumber(capacityName);
1544                        if (Double.isNaN(v) || v < 0)
1545                                v = INFINITE_CAPACITY;
1546                        capacity = (int) v;
1547
1548                        v = edge.getNumber(costName);
1549                        if (Double.isNaN(v))
1550                                v = 1;
1551                        cost = new BigMNumber((int) v);
1552
1553                        String sourceId = edge.getSourceNode().getId();
1554                        String targetId = edge.getTargetNode().getId();
1555                        source = nodes.get(sameDirection ? sourceId : targetId);
1556                        target = nodes.get(sameDirection ? targetId : sourceId);
1557                }
1558
1559                /**
1560                 * Default constructor.
1561                 */
1562                NSArc() {
1563                        cost = new BigMNumber();
1564                }
1565
1566                /**
1567                 * Computes the reduced cost of this arc
1568                 * 
1569                 * @param reducedCost
1570                 *            The result is stored here
1571                 */
1572                void computeReducedCost(BigMNumber reducedCost) {
1573                        reducedCost.set(cost);
1574                        reducedCost.minus(source.potential);
1575                        reducedCost.plus(target.potential);
1576                        if (status == ArcStatus.NONBASIC_UPPER)
1577                                reducedCost.minus();
1578                }
1579
1580                /**
1581                 * Computes the maximum allowed flow change ot this arc.
1582                 * 
1583                 * @param first
1584                 *            One of the endpoints of the arc. Determines the direction
1585                 *            in the cycle.
1586                 * @param flowChange
1587                 *            The result is stored here
1588                 */
1589                void computeAllowedFlowChange(NSNode first, BigMNumber flowChange) {
1590                        if (first == source) {
1591                                // the arc is in the direction of the cycle
1592                                if (capacity == INFINITE_CAPACITY)
1593                                        flowChange.set(0, 1);
1594                                else
1595                                        flowChange.set(capacity - flow);
1596                        } else {
1597                                flowChange.set(flow);
1598                        }
1599                }
1600
1601                /**
1602                 * Changes the flow on this arc.
1603                 * 
1604                 * @param delta
1605                 *            Flow change
1606                 * @param first
1607                 *            One of the endpoints of the arc. Determines the direction
1608                 *            of the cycle
1609                 */
1610                void changeFlow(int delta, NSNode first) {
1611                        if (first == source)
1612                                flow += delta;
1613                        else
1614                                flow -= delta;
1615                        if (animationDelay > 0)
1616                                setUIClass();
1617                }
1618
1619                /**
1620                 * Returns the node on the other side.
1621                 * 
1622                 * @param node
1623                 *            One of the endpoints of this arc
1624                 * @return The opposite node
1625                 */
1626                NSNode getOpposite(NSNode node) {
1627                        if (node == source)
1628                                return target;
1629                        if (node == target)
1630                                return source;
1631                        return null;
1632                }
1633
1634                /**
1635                 * Checks if this arc is artificial.
1636                 * 
1637                 * @return True if the arc is artificial
1638                 */
1639                public boolean isArtificial() {
1640                        return source == root || target == root;
1641                }
1642
1643                /**
1644                 * Returns the id of the edge of the original graph corresponding to
1645                 * this arc.
1646                 * 
1647                 * @return The id of the original edge
1648                 */
1649                String getOriginalId() {
1650                        if (isArtificial())
1651                                return null;
1652                        if (id.startsWith(PREFIX + "REVERSE_"))
1653                                return id.substring(PREFIX.length() + "REVERSE_".length());
1654                        return id;
1655                }
1656
1657                /**
1658                 * Sets the label and the {@code "ui.class"} attribute of the edge (or
1659                 * node if the arc is artificial) of the original graph corresponding to
1660                 * this arc.
1661                 */
1662                void setUIClass() {
1663                        if (isArtificial()) {
1664                                NSNode node = getOpposite(root);
1665                                String uiClass = "trans";
1666                                if (node.supply > 0)
1667                                        uiClass = "supply";
1668                                else if (node.supply < 0)
1669                                        uiClass = "demand";
1670                                uiClass += flow == 0 ? "_balanced" : "_unbalanced";
1671                                Node x = graph.getNode(node.id);
1672                                x.addAttribute("label", target == root ? flow : -flow);
1673                                x.addAttribute("ui.class", uiClass);
1674                        } else {
1675                                String uiClass = "basic";
1676                                if (status == ArcStatus.NONBASIC_LOWER)
1677                                        uiClass = "nonbasic_lower";
1678                                else if (status == ArcStatus.NONBASIC_UPPER)
1679                                        uiClass = "nonbasic_upper";
1680
1681                                Edge e = graph.getEdge(getOriginalId());
1682                                e.addAttribute("label", flow);
1683                                e.addAttribute("ui.class", uiClass);
1684                        }
1685                }
1686
1687                /**
1688                 * Switches arc direction. Use with caution! Used only in
1689                 * {@link NetworkSimplex#changeSupply(NSNode, int)} for artificial basic
1690                 * arcs.
1691                 */
1692                void switchDirection() {
1693                        NSNode tmp = source;
1694                        source = target;
1695                        target = tmp;
1696                        flow = -flow;
1697
1698                        NSNode subtreeRoot = getOpposite(root);
1699                        subtreeRoot.computePotential();
1700                        for (NSNode node = subtreeRoot.thread; node.depth > subtreeRoot.depth; node = node.thread)
1701                                node.computePotential();
1702                }
1703        }
1704
1705        // test and debug
1706
1707        /**
1708         * Prints a table containing informations about the current basic feasible
1709         * solution. Useful for testing and debugging purposes.
1710         * 
1711         * @param ps
1712         *            A stream where the output goes.
1713         */
1714        public void printBFS(PrintStream ps) {
1715                ps.println("=== Nodes ===");
1716                ps.printf("%20s%10s%10s%20s%20s%10s%n", "id", "supply", "potential",
1717                                "parent", "thread", "depth");
1718                ps.printf("%20s%10d%10s%20s%20s%10d%n", root.id, root.supply,
1719                                root.potential, "-", root.thread.id, root.depth);
1720                for (NSNode node : nodes.values())
1721                        ps.printf("%20s%10d%10s%20s%20s%10d%n", node.id, node.supply,
1722                                        node.potential, node.parent.id, node.thread.id, node.depth);
1723                ps.println();
1724
1725                ps.println("=== Arcs ===");
1726                ps.printf("%20s%10s%10s%10s%10s%20s%n", "id", "capacity", "cost",
1727                                "flow", "r. cost", "status");
1728                for (NSArc a : arcs.values()) {
1729                        a.computeReducedCost(work1);
1730                        ps.printf("%20s%10s%10s%10s%10s%20s%n", a.id,
1731                                        a.capacity == INFINITE_CAPACITY ? "Inf" : a.capacity,
1732                                        a.cost, a.flow, work1, a.status);
1733                }
1734
1735                for (NSNode node : nodes.values()) {
1736                        NSArc a = node.artificialArc;
1737                        a.computeReducedCost(work1);
1738                        ps.printf("%20s%10s%10s%10s%10s%20s%n", a.id,
1739                                        a.capacity == INFINITE_CAPACITY ? "Inf" : a.capacity,
1740                                        a.cost, a.flow, work1, a.status);
1741                }
1742
1743                ps.println();
1744                ps.printf("=== Objective value %s. Solution status %s ===%n%n",
1745                                objectiveValue, solutionStatus);
1746        }
1747}