diff options
| -rw-r--r-- | impl/Complete.hpp | 10 | ||||
| -rw-r--r-- | impl/EquationSystem.g | 18 | ||||
| -rwxr-xr-x | impl/MCFClass.h | 2438 | ||||
| -rw-r--r-- | impl/MCFSimplex.cpp | 4198 | ||||
| -rw-r--r-- | impl/MCFSimplex.h | 1086 | ||||
| -rw-r--r-- | impl/MCFSimplex.o | bin | 0 -> 205412 bytes | |||
| -rw-r--r-- | impl/Makefile | 17 | ||||
| -rwxr-xr-x | impl/OPTUtils.h | 475 | ||||
| -rw-r--r-- | impl/Operator.hpp | 82 | ||||
| -rw-r--r-- | impl/main.cpp | 73 |
10 files changed, 8383 insertions, 14 deletions
diff --git a/impl/Complete.hpp b/impl/Complete.hpp index e3ec15a..664d71f 100644 --- a/impl/Complete.hpp +++ b/impl/Complete.hpp @@ -99,6 +99,16 @@ struct Complete { template<typename Z> friend std::ostream& operator<<(std::ostream&, const Complete<Z>&); + template<typename S> + S as() const { + if (_infinity) { + if (_value > 0) + return infinity<S>(); + return -infinity<S>(); + } + return (S) _value; + } + private: T _value; bool _infinity; diff --git a/impl/EquationSystem.g b/impl/EquationSystem.g index 4884ca9..3ef9009 100644 --- a/impl/EquationSystem.g +++ b/impl/EquationSystem.g @@ -14,6 +14,7 @@ tokens { GUARD = 'guard' ; GREATER_EQUAL = '>=' ; QUESTION_MARK = '?' ; + MCF = 'MCF' ; MAXIMUM = 'max' ; MINIMUM = 'min' ; NEWLINE = '\n' ; @@ -26,18 +27,23 @@ maxExpr : MAXIMUM^ '('! minExpr ( ','! minExpr )* ')'! | minExpr ; minExpr : MINIMUM^ '('! maxExpr ( ','! maxExpr )* ')'! | expr ; expr : '(' expr GREATER_EQUAL expr QUESTION_MARK expr ')' -> ^(GUARD expr expr expr) - | GUARD^ '('! maxExpr ','! maxExpr ','! maxExpr ')'! - | 'add'^ '('! maxExpr ( ','! maxExpr )* ')'! - | 'sub'^ '('! maxExpr ( ','! maxExpr )* ')'! - | 'mult'^ '('! maxExpr ( ','! maxExpr )* ')'! - | term ( (PLUS | MULT | SUB | COMMA)^ expr )* ; + | GUARD^ '('! maxExpr ','! maxExpr ','! maxExpr ')'! + | MCF^ '<'! supplies ','! arcs '>'! '('! maxExpr (','! maxExpr)* ')'! + | 'add'^ '('! maxExpr ( ','! maxExpr )* ')'! + | 'sub'^ '('! maxExpr ( ','! maxExpr )* ')'! + | 'mult'^ '('! maxExpr ( ','! maxExpr )* ')'! + | term ( (PLUS | MULT | SUB | COMMA)^ expr )* ; + +supplies : '['^ NUMBER (','! NUMBER)* ']'! ; +arc : NUMBER ':'^ NUMBER ; +arcs : '['^ arc (','! arc)* ']'! ; term : NUMBER | VARIABLE | '-'^ term ; -NUMBER : (DIGIT)+ ; +NUMBER : '-'? (DIGIT)+ ; VARIABLE: (LETTER) (LETTER | DIGIT | '-' | '\[' | ']' )* ; WHITESPACE : ( '\t' | ' ' | '\u000C' )+ { diff --git a/impl/MCFClass.h b/impl/MCFClass.h new file mode 100755 index 0000000..f10cb5a --- /dev/null +++ b/impl/MCFClass.h @@ -0,0 +1,2438 @@ +/*--------------------------------------------------------------------------*/ +/*-------------------------- File MCFClass.h -------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @file + * Header file for the abstract base class MCFClass, which defines a standard + * interface for (linear or convex quadratic separable) Min Cost Flow Problem + * solvers, to be implemented as derived classes. + * + * \version 3.01 + * + * \date 30 - 09 - 2011 + * + * \author Alessandro Bertolini \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * \author Antonio Frangioni \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * \author Claudio Gentile \n + * Istituto di Analisi di Sistemi e Informatica \n + * Consiglio Nazionale delle Ricerche \n + * + * Copyright © 1996 - 2011 by Antonio Frangioni, Claudio Gentile + */ + +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*----------------------------- DEFINITIONS --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef __MCFClass + #define __MCFClass /* self-identification: #endif at the end of the file */ + +/*--------------------------------------------------------------------------*/ +/*--------------------------------- MACROS ---------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup MCFCLASS_MACROS Compile-time switches in MCFClass.h + These macros control some important details of the class interface. + Although using macros for activating features of the interface is not + very C++, switching off some unused features may allow some + implementation to be more efficient in running time or memory. + @{ */ + +/*-------------------------------- USENAME0 --------------------------------*/ + +#define USENAME0 1 + +/**< Decides if 0 or 1 is the "name" of the first node. + If USENAME0 == 1, (warning: it has to be *exactly* 1), then the node + names go from 0 to n - 1, otherwise from 1 to n. Note that this does not + affect the position of the deficit in the deficit vectors, i.e., the + deficit of the i-th node - be its "name" `i' or `i - 1' - is always in + the i-th position of the vector. */ + +/*@} end( group( MCFCLASS_MACROS ) ) */ +/*--------------------------------------------------------------------------*/ +/*------------------------------ INCLUDES ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#include "OPTUtils.h" + +/* OPTtypes.h defines standard interfaces for timing and random routines, as + well as the namespace OPTtypes_di_unipi_it and the macro + OPT_USE_NAMESPACES, useful for switching off all namespaces in one blow + for those strange cases where they create problems. */ + +#include <iomanip> +#include <sstream> +#include <limits> +#include <cassert> + +// QUICK HACK +#define throw(x) assert(false) + +/*--------------------------------------------------------------------------*/ +/*------------------------- NAMESPACE and USING ----------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) +namespace MCFClass_di_unipi_it +{ + /** @namespace MCFClass_di_unipi_it + The namespace MCFClass_di_unipi_it is defined to hold the MCFClass + class and all the relative stuff. It comprises the namespace + OPTtypes_di_unipi_it. */ + + using namespace OPTtypes_di_unipi_it; +#endif + +/*@} end( group( MCFCLASS_CONSTANTS ) ) */ +/*--------------------------------------------------------------------------*/ +/*-------------------------- CLASS MCFClass --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup MCFCLASS_CLASSES Classes in MCFClass.h + @{ */ + +/*--------------------------------------------------------------------------*/ +/*--------------------------- GENERAL NOTES --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** This abstract base class defines a standard interface for (linear or + convex quadartic separable) Min Cost Flow (MCF) problem solvers. + + The data of the problem consist of a (directed) graph G = ( N , A ) with + n = |N| nodes and m = |A| (directed) arcs. Each node `i' has a deficit + b[ i ], i.e., the amount of flow that is produced/consumed by the node: + source nodes (which produce flow) have negative deficits and sink nodes + (which consume flow) have positive deficits. Each arc `(i, j)' has an + upper capacity U[ i , j ], a linear cost coefficient C[ i , j ] and a + (non negative) quadratic cost coefficient Q[ i , j ]. Flow variables + X[ i , j ] represents the amount of flow to be sent on arc (i, j). + Parallel arcs, i.e., multiple copies of the same arc `(i, j)' (with + possibily different costs and/or capacities) are in general allowed. + The formulation of the problem is therefore: + \f[ + \min \sum_{ (i, j) \in A } C[ i , j ] X[ i, j ] + + Q[ i , j ] X[ i, j ]^2 / 2 + \f] + \f[ + (1) \sum_{ (j, i) \in A } X[ j , i ] - + \sum_{ (i, j) \in A } X[ i , j ] = b[ i ] + \hspace{1cm} i \in N + \f] + \f[ + (2) 0 \leq X[ i , j ] \leq U[ i , j ] + \hspace{1cm} (i, j) \in A + \f] + The n equations (1) are the flow conservation constraints and the 2m + inequalities (2) are the flow nonnegativity and capacity constraints. + At least one of the flow conservation constraints is redundant, as the + demands must be balanced (\f$\sum_{ i \in N } b[ i ] = 0\f$); indeed, + exactly n - ConnectedComponents( G ) flow conservation constraints are + redundant, as demands must be balanced in each connected component of G. + Let us denote by QA and LA the disjoint subsets of A containing, + respectively, "quadratic" arcs (with Q[ i , j ] > 0) and "linear" arcs + (with Q[ i , j ] = 0); the (MCF) problem is linear if QA is empty, and + nonlinear (convex quadratic) if QA is nonempty. + + The dual of the problem is: + \f[ + \max \sum_{ i \in N } Pi[ i ] b[ i ] - + \sum_{ (i, j) \in A } W[ i , j ] U[ i , j ] - + \sum_{ (i, j) \in AQ } V[ i , j ]^2 / ( 2 * Q[ i , j ] ) + \f] + \f[ + (3.a) C[ i , j ] - Pi[ j ] + Pi[ i ] + W[ i , j ] - Z[ i , j ] = 0 + \hspace{1cm} (i, j) \in AL + \f] + \f[ + (3.b) C[ i , j ] - Pi[ j ] + Pi[ i ] + W[ i , j ] - Z[ i , j ] = + V[ i , j ] + \hspace{1cm} (i, j) \in AQ + \f] + \f[ + (4.a) W[ i , j ] \geq 0 \hspace{1cm} (i, j) \in A + \f] + \f[ + (4.b) Z[ i , j ] \geq 0 \hspace{1cm} (i, j) \in A + \f] + + Pi[] is said the vector of node potentials for the problem, W[] are + bound variables and Z[] are slack variables. Given Pi[], the quantities + \f[ + RC[ i , j ] = C[ i , j ] + Q[ i , j ] * X[ i , j ] - Pi[ j ] + Pi[ i ] + \f] + are said the "reduced costs" of arcs. + + A primal and dual feasible solution pair is optimal if and only if the + complementary slackness conditions + \f[ + RC[ i , j ] > 0 \Rightarrow X[ i , j ] = 0 + \f] + \f[ + RC[ i , j ] < 0 \Rightarrow X[ i , j ] = U[ i , j ] + \f] + are satisfied for all arcs (i, j) of A. + + The MCFClass class provides an interface with methods for managing and + solving problems of this kind. Actually, the class can also be used as + an interface for more general NonLinear MCF problems, where the cost + function either nonseparable ( C( X ) ) or arc-separable + ( \f$\sum_{ (i, j) \in A } C_{i,j}( X[ i, j ] )\f$ ). However, solvers + for NonLinear MCF problems are typically objective-function-specific, + and there is no standard way for inputting a nonlinear function different + from a separable convex quadratic one, so only the simplest form is dealt + with in the interface, leaving more complex NonLinear parts to the + interface of derived classes. */ + +class MCFClass { + +/*--------------------------------------------------------------------------*/ +/*----------------------- PUBLIC PART OF THE CLASS -------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-- --*/ +/*-- The following methods and data are the actual interface of the --*/ +/*-- class: the standard user should use these methods and data only. --*/ +/*-- --*/ +/*--------------------------------------------------------------------------*/ + + public: + +/*--------------------------------------------------------------------------*/ +/*---------------------------- PUBLIC TYPES --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Public types + The MCFClass defines four main public types: + + - Index, the type of arc and node indices; + + - FNumber, the type of flow variables, arc capacities, and node deficits; + + - CNumber, the type of flow costs, node potentials, and arc reduced costs; + + - FONumber, the type of objective function value. + + By re-defining the types in this section, most MCFSolver should be made + to work with any reasonable choice of data type (= one that is capable of + properly representing the data of the instances to be solved). This may + be relevant due to an important property of MCF problems: *if all arc + capacities and node deficits are integer, then there exists an integral + optimal primal solution*, and *if all arc costs are integer, then there + exists an integral optimal dual solution*. Even more importantly, *many + solution algorithms will in fact produce an integral primal/dual + solution for free*, because *every primal/dual solution they generate + during the solution process is naturally integral*. Therefore, one can + use integer data types to represent everything connected with flows and/or + costs if the corresponding data is integer in all instances one needs to + solve. This directly translates in significant memory savings and/or speed + improvements. + + *It is the user's responsibility to ensure that these types are set to + reasonable values*. So, the experienced user may want to experiment with + setting this data properly if memory footprint and/or speed is a primary + concern. Note, however, that *not all solution algorithms will happily + accept integer data*; one example are Interior-Point approaches, which + require both flow and cost variables to be continuous (float). So, the + viability of setting integer data (as well as its impact on performances) + is strictly related to the specific kind of algorithm used. Since these + types are common to all derived classes, they have to be set taking into + account the needs of all the solvers that are going to be used, and + adapting to the "worst case"; of course, FNumber == CNumber == double is + going to always be an acceptable "worst case" setting. MCFClass may in a + future be defined as a template class, with these as template parameters, + but this is currently deemed overkill and avoided. + + Finally, note that the above integrality property only holds for *linear* + MCF problems. If any arc has a nonzero quadratic cost coefficient, optimal + flows and potentials may be fractional even if all the data of the problem + (comprised quadratic cost coefficients) is integer. Hence, for *quadratic* + MCF solvers, a setting like FNumber == CNumber == double is actually + *mandatory*, for any reasonable algorithm will typically misbehave + otherwise. + @{ */ + +/*--------------------------------------------------------------------------*/ + + typedef unsigned int Index; ///< index of a node or arc ( >= 0 ) + typedef Index *Index_Set; ///< set (array) of indices + typedef const Index cIndex; ///< a read-only index + typedef cIndex *cIndex_Set; ///< read-only index array + +/*--------------------------------------------------------------------------*/ + + typedef int SIndex; ///< index of a node or arc + typedef SIndex *SIndex_Set; ///< set (array) of indices + typedef const SIndex cSIndex; ///< a read-only index + typedef cSIndex *cSIndex_Set; ///< read-only index array + +/*--------------------------------------------------------------------------*/ + + typedef double FNumber; ///< type of arc flow + typedef FNumber *FRow; ///< vector of flows + typedef const FNumber cFNumber; ///< a read-only flow + typedef cFNumber *cFRow; ///< read-only flow array + +/*--------------------------------------------------------------------------*/ + + typedef double CNumber; ///< type of arc flow cost + typedef CNumber *CRow; ///< vector of costs + typedef const CNumber cCNumber; ///< a read-only cost + typedef cCNumber *cCRow; ///< read-only cost array + +/*--------------------------------------------------------------------------*/ + + typedef double FONumber; + /**< type of the objective function: has to hold sums of products of + FNumber(s) by CNumber(s) */ + + typedef const FONumber cFONumber; ///< a read-only o.f. value + +/*--------------------------------------------------------------------------*/ +/** Very small class to simplify extracting the "+ infinity" value for a + basic type (FNumber, CNumber, Index); just use Inf<type>(). */ + + template <typename T> + class Inf { + public: + Inf() {} + operator T() { return( std::numeric_limits<T>::max() ); } + }; + +/*--------------------------------------------------------------------------*/ +/** Very small class to simplify extracting the "machine epsilon" for a + basic type (FNumber, CNumber); just use Eps<type>(). */ + + template <typename T> + class Eps { + public: + Eps() {} + operator T() { return( std::numeric_limits<T>::epsilon() ); } + }; + +/*--------------------------------------------------------------------------*/ +/** Small class for exceptions. Derives from std::exception implementing the + virtual method what() -- and since what is virtual, remember to always + catch it by reference (catch exception &e) if you want the thing to work. + MCFException class are thought to be of the "fatal" type, i.e., problems + for which no solutions exists apart from aborting the program. Other kinds + of exceptions can be properly handled by defining derived classes with + more information. */ + + class MCFException : public exception { + public: + MCFException( const char *const msg = 0 ) { errmsg = msg; } + + // wow, such a hack +#undef throw + const char* what( void ) const throw () { return( errmsg ); } +#define throw(x) assert(false) + private: + const char *errmsg; + }; + +/*--------------------------------------------------------------------------*/ +/** Public enum describing the possible parameters of the MCF solver, to be + used with the methods SetPar() and GetPar(). */ + + enum MCFParam { kMaxTime = 0 , ///< max time + kMaxIter , ///< max number of iteration + kEpsFlw , ///< tolerance for flows + kEpsDfct , ///< tolerance for deficits + kEpsCst , ///< tolerance for costs + kReopt , ///< whether or not to reoptimize + kLastParam /**< dummy parameter: this is used to + allow derived classes to "extend" + the set of parameters. */ + }; + +/*--------------------------------------------------------------------------*/ +/** Public enum describing the possible status of the MCF solver. */ + + enum MCFStatus { kUnSolved = -1 , ///< no solution available + kOK = 0 , ///< optimal solution found + + kStopped , ///< optimization stopped + kUnfeasible , ///< problem is unfeasible + kUnbounded , ///< problem is unbounded + kError ///< error in the solver + }; + +/*--------------------------------------------------------------------------*/ +/** Public enum describing the possible reoptimization status of the MCF + solver. */ + + enum MCFAnswer { kNo = 0 , ///< no + kYes ///< yes + }; + +/*--------------------------------------------------------------------------*/ +/** Public enum describing the possible file formats in WriteMCF(). */ + + enum MCFFlFrmt { kDimacs = 0 , ///< DIMACS file format for MCF + kQDimacs , ///< quadratic DIMACS file format for MCF + kMPS , ///< MPS file format for LP + kFWMPS ///< "Fixed Width" MPS format + }; + +/*--------------------------------------------------------------------------*/ +/** Base class for representing the internal state of the MCF algorithm. */ + + class MCFState { + public: + MCFState( void ) {}; + virtual ~MCFState() {}; + }; + + typedef MCFState *MCFStatePtr; ///< pointer to a MCFState + +/*@} -----------------------------------------------------------------------*/ +/*--------------------------- PUBLIC METHODS -------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*---------------------------- CONSTRUCTOR ---------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Constructors + @{ */ + + MCFClass( cIndex nmx = 0 , cIndex mmx = 0 ) + { + nmax = nmx; + mmax = mmx; + n = m = 0; + + status = kUnSolved; + Senstv = true; + + EpsFlw = Eps<FNumber>() * 100; + EpsCst = Eps<CNumber>() * 100; + EpsDfct = EpsFlw * ( nmax ? nmax : 100 ); + + MaxTime = 0; + MaxIter = 0; + + MCFt = NULL; + } + +/**< Constructor of the class. + + nmx and mmx, if provided, are taken to be respectively the maximum number + of nodes and arcs in the network. If nonzero values are passed, memory + allocation can be anticipated in the constructor, which is sometimes + desirable. The maximum values are stored in the protected fields nmax and + mmax, and can be changed with LoadNet() [see below]; however, changing + them typically requires memory allocation/deallocation, which is + sometimes undesirable outside the constructor. + + After that an object has been constructed, no problem is loaded; this has + to be done with LoadNet() [see below]. Thus, it is an error to invoke any + method which requires the presence of a problem (typicall all except those + in the initializations part). The base class provides two protected fields + n and m for the current number of nodes and arcs, respectively, that are + set to 0 in the constructor precisely to indicate that no instance is + currently loaded. */ + +/*@} -----------------------------------------------------------------------*/ +/*-------------------------- OTHER INITIALIZATIONS -------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Other initializations + @{ */ + + virtual void LoadNet( cIndex nmx = 0 , cIndex mmx = 0 , cIndex pn = 0 , + cIndex pm = 0 , cFRow pU = NULL , cCRow pC = NULL , + cFRow pDfct = NULL , cIndex_Set pSn = NULL , + cIndex_Set pEn = NULL ) = 0; + +/**< Inputs a new network. + + The parameters nmx and mmx are the new max number of nodes and arcs, + possibly overriding those set in the constructor [see above], altough at + the likely cost of memory allocation and deallocation. Passing nmx == + mmx == 0 is intended as a signal to the solver to deallocate everything + and wait for new orders; in this case, all the other parameters are ignored. + + Otherwise, in principle all the other parameters have to be provided. + Actually, some of them may not be needed for special classes of MCF + problems (e.g., costs in a MaxFlow problem, or start/end nodes in a + problem defined over a graph with fixed topology, such as a complete + graph). Also, passing NULL is allowed to set default values. + + The meaning of the parameters is the following: + + - pn is the current number of nodes of the network (<= nmax). + - pm is the number of arcs of the network (<= mmax). + + - pU is the m-vector of the arc upper capacities; capacities must be + nonnegative, but can in principle be infinite (== F_INF); passing + pU == NULL means that all capacities are infinite; + + - pC is the m-vector of the arc costs; costs must be finite (< C_INF); + passing pC == NULL means that all costs must be 0. + + - pDfct is the n-vector of the node deficits; source nodes have negative + deficits and sink nodes have positive deficits; passing pDfct == + NULL means that all deficits must be 0 (a circulation problem); + + - pSn is the m-vector of the arc starting nodes; pSn == NULL is in + principle not allowed, unless the topology of the graph is fixed; + - pEn is the m-vector of the arc ending nodes; same comments as for pSn. + + Note that node "names" in the arrays pSn and pEn must go from 1 to pn if + the macro USANAME0 [see above] is set to 0, while they must go from 0 to + pn - 1 if USANAME0 is set to 1. In both cases, however, the deficit of the + first node is read from the first (0-th) position of pDfct, that is if + USANAME0 == 0 then the deficit of the node with name `i' is read from + pDfct[ i - 1 ]. + + The data passed to LoadNet() can be used to specify that the arc `i' must + not "exist" in the problem. This is done by passing pC[ i ] == C_INF; + solvers which don't read costs are forced to read them in order to check + this, unless they provide alternative solver-specific ways to accomplish + the same tasks. These arcs are "closed", as for the effect of CloseArc() + [see below]. "invalid" costs (== C_INF) are set to 0 in order to being + subsequently capable of "opening" them back with OpenArc() [see below]. + The way in which these non-existent arcs are phisically dealt with is + solver-specific; in some solvers, for instance, this could be obtained by + simply putting their capacity to zero. Details about these issues should + be found in the interface of derived classes. + + Note that the quadratic part of the objective function, if any, is not + dealt with in LoadNet(); it can only be separately provided with + ChgQCoef() [see below]. By default, the problem is linear, i.e., all + coefficients of the second-order terms in the objective function are + assumed to be zero. */ + +/*--------------------------------------------------------------------------*/ + + virtual inline void LoadDMX( istream &DMXs , bool IsQuad = false ); + +/**< Read a MCF instance in DIMACS standard format from the istream. The + format is the following. The first line must be + + p min <number of nodes> <number of arcs> + + Then the node definition lines must be found, in the form + + n <node number> <node supply> + + Not all nodes need have a node definition line; these are given zero + supply, i.e., they are transhipment nodes (supplies are the inverse of + deficits, i.e., a node with positive supply is a source node). Finally, + the arc definition lines must be found, in the form + + a <start node> <end node> <lower bound> <upper bound> <flow cost> + + There must be exactly <number of arcs> arc definition lines in the file. + + This method is *not* pure virtual because an implementation is provided by + the base class, using the LoadNet() method (which *is* pure virtual). + However, the method *is* virtual to allow derived classes to implement + more efficient versions, should they have any reason to do so. + + \note Actually, the file format accepted by LoadDMX (at least in the + base class implementation) is more general than the DIMACS standard + format, in that it is allowed to mix node and arc definitions in + any order, while the DIMACS file requires all node information to + appear before all arc information. + + \note Other than for the above, this method is assumed to allow for + *quadratic* Dimacs files, encoding for convex quadratic separable + Min Cost Flow instances. This is a simple extension where each arc + descriptor has a sixth field, <quadratic cost>. The provided + istream is assumed to be quadratic Dimacs file if IsQuad is true, + and a regular linear Dimacs file otherwise. */ + +/*--------------------------------------------------------------------------*/ + + virtual void PreProcess( void ) {} + +/**< Extract a smaller/easier equivalent MCF problem. The data of the instance + is changed and the easier one is solved instead of the original one. In the + MCF case, preprocessing may involve reducing bounds, identifying + disconnected components of the graph etc. However, proprocessing is + solver-specific. + + This method can be implemented by derived classes in their solver-specific + way. Preprocessing may reveal unboundedness or unfeasibility of the + problem; if that happens, PreProcess() should properly set the `status' + field, that can then be read with MCFGetStatus() [see below]. + + Note that preprocessing may destroy all the solution information. Also, it + may be allowed to change the data of the problem, such as costs/capacities + of the arcs. + + A valid preprocessing is doing nothing, and that's what the default + implementation of this method (that is *not* pure virtual) does. */ + +/*--------------------------------------------------------------------------*/ + + virtual inline void SetPar( int par , int val ); + +/**< Set integer parameters of the algorithm. + + @param par is the parameter to be set; the enum MCFParam can be used, but + 'par' is an int (every enum is an int) so that the method can + be extended by derived classes for the setting of their + parameters + + @param value is the value to assign to the parameter. + + The base class implementation handles these parameters: + + - kMaxIter: the max number of iterations in which the MCF Solver can find + an optimal solution (default 0, which means no limit) + + - kReopt: tells the solver if it has to reoptimize. The implementation in + the base class sets a flag, the protected \c bool field \c + Senstv; if true (default) this field instructs the MCF solver to + to try to exploit the information about the latest optimal + solution to speedup the optimization of the current problem, + while if the field is false the MCF solver should restart the + optimization "from scratch" discarding any previous information. + Usually reoptimization speeds up the computation considerably, + but this is not always true, especially if the data of the + problem changes a lot.*/ + +/*--------------------------------------------------------------------------*/ + + virtual inline void SetPar( int par , double val ); + +/**< Set float parameters of the algorithm. + + @param par is the parameter to be set; the enum MCFParam can be used, but + 'par' is an int (every enum is an int) so that the method can + be extended by derived classes for the setting of their + parameters + + @param value is the value to assign to the parameter. + + The base class implementation handles these parameters: + + - kEpsFlw: sets the tolerance for controlling if the flow on an arc is zero + to val. This also sets the tolerance for controlling if a node + deficit is zero (see kEpsDfct) to val * < max number of nodes >; + this value should be safe for graphs in which any node has less + than < max number of nodes > adjacent nodes, i.e., for all graphs + but for very dense ones with "parallel arcs" + + - kEpsDfct: sets the tolerance for controlling if a node deficit is zero to + val, in case a better value than that autmatically set by + kEpsFlw (see above) is available (e.g., val * k would be good + if no node has more than k neighbours) + + - kEpsCst: sets the tolerance for controlling if the reduced cost of an arc + is zero to val. A feasible solution satisfying eps-complementary + slackness, i.e., such that + \f[ + RC[ i , j ] < - eps \Rightarrow X[ i , j ] = U[ ij ] + \f] + and + \f[ + RC[ i , j ] > eps \Rightarrow X[ i , j ] == 0 , + \f] + is known to be ( eps * n )-optimal. + + - kMaxTime: sets the max time (in seconds) in which the MCF Solver can find + an optimal solution (default 0, which means no limit). */ + +/*--------------------------------------------------------------------------*/ + + virtual inline void GetPar( int par , int &val ); + +/**< This method returns one of the integer parameter of the algorithm. + + @param par is the parameter to return [see SetPar( int ) for comments]; + + @param val upon return, it will contain the value of the parameter. + + The base class implementation handles the parameters kMaxIter and kReopt. + */ + +/*--------------------------------------------------------------------------*/ + + virtual inline void GetPar( int par , double &val ); + +/**< This method returns one of the integer parameter of the algorithm. + + @param par is the parameter to return [see SetPar( double ) for comments]; + + @param val upon return, it will contain the value of the parameter. + + The base class implementation handles the parameters kEpsFlw, kEpsDfct, + kEpsCst, and kMaxTime. */ + +/*--------------------------------------------------------------------------*/ + + virtual void SetMCFTime( bool TimeIt = true ) + { + if( TimeIt ) + if( MCFt ) + MCFt->ReSet(); + else + MCFt = new OPTtimers(); + else { + delete MCFt; + MCFt = NULL; + } + } + +/**< Allocate an OPTtimers object [see OPTtypes.h] to be used for timing the + methods of the class. The time can be read with TimeMCF() [see below]. By + default, or if SetMCFTime( false ) is called, no timing is done. Note that, + since all the relevant methods ot the class are pure virtual, MCFClass can + only manage the OPTtimers object, but it is due to derived classes to + actually implement the timing. + + @note time accumulates over the calls: calling SetMCFTime(), however, + resets the counters, allowing to time specific groups of calls. + + @note of course, setting kMaxTime [see SetPar() above] to any nonzero + value has no effect unless SetMCFTime( true ) has been called. */ + +/*@} -----------------------------------------------------------------------*/ +/*---------------------- METHODS FOR SOLVING THE PROBLEM -------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Solving the problem + @{ */ + + virtual void SolveMCF( void ) = 0; + +/**< Solver of the Min Cost Flow Problem. Attempts to solve the MCF instance + currently loaded in the object. */ + +/*--------------------------------------------------------------------------*/ + + inline int MCFGetStatus( void ) + { + return( status ); + } + +/**< Returns an int describing the current status of the MCF solver. Possible + return values are: + + - kUnSolved SolveMCF() has not been called yet, or the data of the + problem has been changed since the last call; + + - kOK optimization has been carried out succesfully; + + - kStopped optimization have been stopped before that the stopping + conditions of the solver applied, e.g. because of the + maximum allowed number of "iterations" [see SetPar( int )] + or the maximum allowed time [see SetPar( double )] has been + reached; this is not necessarily an error, as it might just + be required to re-call SolveMCF() giving it more "resources" + in order to solve the problem; + + - kUnfeasible if the current MCF instance is (primal) unfeasible; + + - kUnbounded if the current MCF instance is (primal) unbounded (this can + only happen if the solver actually allows F_INF capacities, + which is nonstandard in the interface); + + - kError if there was an error during the optimization; this + typically indicates that computation cannot be resumed, + although solver-dependent ways of dealing with + solver-dependent errors may exist. + + MCFClass has a protected \c int \c member \c status \c that can be used by + derived classes to hold status information and that is returned by the + standard implementation of this method. Note that \c status \c is an + \c int \c and not an \c enum \c, and that an \c int \c is returned by this + method, in order to allow the derived classes to extend the set of return + values if they need to do so. */ + +/*@} -----------------------------------------------------------------------*/ +/*---------------------- METHODS FOR READING RESULTS -----------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Reading flow solution + @{ */ + + virtual void MCFGetX( FRow F , Index_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Write the optimal flow solution in the vector F[]. If nms == NULL, F[] + will be in "dense" format, i.e., the flow relative to arc `i' + (i in 0 .. m - 1) is written in F[ i ]. If nms != NULL, F[] will be in + "sparse" format, i.e., the indices of the nonzero elements in the flow + solution are written in nms (that is then Inf<Index>()-terminated) and + the flow value of arc nms[ i ] is written in F[ i ]. Note that nms is + *not* guaranteed to be ordered. Also, note that, unlike MCFGetRC() and + MCFGetPi() [see below], nms is an *output* of the method. + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cFRow MCFGetX( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal data structure containing the + flow solution in "dense" format. Since this may *not always be available*, + depending on the implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual bool HaveNewX( void ) + { + return( false ); + } + +/**< Return true if a different (approximately) optimal primal solution is + available. If the method returns true, then any subsequent call to (any + form of) MCFGetX() will return a different primal solution w.r.t. the one + that was being returned *before* the call to HaveNewX(). This solution need + not be optimal (although, ideally, it has to be "good); this can be checked + by comparing its objective function value, that will be returned by a call + to MCFGetFO() [see below]. + + Any subsequent call of HaveNewX() that returns true produces a new + solution, until the first that returns false; from then on, no new + solutions will be generated until something changes in the problem's + data. + + Note that a default implementation of HaveNewX() is provided which is + good for those solvers that only produce one optimal primal solution. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading potentials + @{ */ + + virtual void MCFGetPi( CRow P , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Writes the optimal node potentials in the vector P[]. If nms == NULL, + the node potential of node `i' (i in 0 .. n - 1) is written in P[ i ] + (note that here node names always start from zero, regardless to the value + of USENAME0). If nms != NULL, it must point to a vector of indices in + 0 .. n - 1 (ordered in increasing sense and Inf<Index>()-terminated), and + the node potential of nms[ i ] is written in P[ i ]. Note that, unlike + MCFGetX() above, nms is an *input* of the method. + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the nodes `i' with strt <= i < min( MCFn() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to nodes which are *both* in nms[] and whose index is + in the correct range are returned. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cCRow MCFGetPi( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal data structure containing the + node potentials. Since this may *not always be available*, depending on + the implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual bool HaveNewPi( void ) + { + return( false ); + } + +/**< Return true if a different (approximately) optimal dual solution is + available. If the method returns true, then any subsequent call to (any + form of) MCFGetPi() will return a different dual solution, and MCFGetRC() + [see below] will return the corresponding reduced costs. The new solution + need not be optimal (although, ideally, it has to be "good); this can be + checked by comparing its objective function value, that will be returned + by a call to MCFGetDFO() [see below]. + + Any subsequent call of HaveNewPi() that returns true produces a new + solution, until the first that returns false; from then on, no new + solutions will be generated until something changes in the problem's + data. + + Note that a default implementation of HaveNewPi() is provided which is + good for those solvers that only produce one optimal dual solution. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading reduced costs + @{ */ + + virtual void MCFGetRC( CRow CR , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Write the reduced costs corresponding to the current dual solution in + RC[]. If nms == NULL, the reduced cost of arc `i' (i in 0 .. m - 1) is + written in RC[ i ]; if nms != NULL, it must point to a vector of indices + in 0 .. m - 1 (ordered in increasing sense and Inf<Index>()-terminated), + and the reduced cost of arc nms[ i ] is written in RC[ i ]. Note that, + unlike MCFGetX() above, nms is an *input* of the method. + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to arcs which are *both* in nms[] and whose index is + in the correct range are returned. + + @note the output of MCFGetRC() will change after any call to HaveNewPi() + [see above] which returns true. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cCRow MCFGetRC( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal data structure containing the + reduced costs. Since this may *not always be available*, depending on the + implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. + + @note the output of MCFGetRC() will change after any call to HaveNewPi() + [see above] which returns true. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual CNumber MCFGetRC( cIndex i ) = 0; + +/**< Return the reduced cost of the i-th arc. This information should be + cheapily available in most implementations. + + @note the output of MCFGetRC() will change after any call to HaveNewPi() + [see above] which returns true. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading the objective function value + @{ */ + + virtual FONumber MCFGetFO( void ) = 0; + +/**< Return the objective function value of the primal solution currently + returned by MCFGetX(). + + If MCFGetStatus() == kOK, this is guaranteed to be the optimal objective + function value of the problem (to within the optimality tolerances), but + only prior to any call to HaveNewX() that returns true. MCFGetFO() + typically returns Inf<FONumber>() if MCFGetStatus() == kUnfeasible and + - Inf<FONumber>() if MCFGetStatus() == kUnbounded. If MCFGetStatus() == + kStopped and MCFGetFO() returns a finite value, it must be an upper bound + on the optimal objective function value (typically, the objective function + value of one primal feasible solution). */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual FONumber MCFGetDFO( void ) + { + switch( MCFGetStatus() ) { + case( kUnSolved ): + case( kStopped ): + case( kError ): return( -Inf<FONumber>() ); + default: return( MCFGetFO() ); + } + } + +/**< Return the objective function value of the dual solution currently + returned by MCFGetPi() / MCFGetRC(). This value (possibly) changes after + any call to HaveNewPi() that returns true. The relations between + MCFGetStatus() and MCFGetDFO() are analogous to these of MCFGetFO(), + except that a finite value corresponding to kStopped must be a lower + bound on the optimal objective function value (typically, the objective + function value one dual feasible solution). + + A default implementation is provided for MCFGetDFO(), which is good for + MCF solvers where the primal and dual optimal solution values always + are identical (except if the problem is unfeasible/unbounded). */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Getting unfeasibility certificate + @{ */ + + virtual FNumber MCFGetUnfCut( Index_Set Cut ) + { + *Cut = Inf<Index>(); + return( 0 ); + } + +/**< Return an unfeasibility certificate. In an unfeasible MCF problem, + unfeasibility can always be reduced to the existence of a cut (subset of + nodes of the graph) such as either: + + - the inverse of the deficit of the cut (the sum of all the deficits of the + nodes in the cut) is larger than the forward capacity of the cut (sum of + the capacities of forward arcs in the cut); that is, the nodes in the + cut globally produce more flow than can be routed to sinks outside the + cut; + + - the deficit of the cut is larger than the backward capacity of the cut + (sum of the capacities of backward arcs in the cut); that is, the nodes + in the cut globally require more flow than can be routed to them from + sources outside the cut. + + When detecting unfeasibility, MCF solvers are typically capable to provide + one such cut. This information can be useful - typically, the only way to + make the problem feasible is to increase the capacity of at least one of + the forward/backward arcs of the cut -, and this method is provided for + getting it. It can be called only if MCFGetStatus() == kUnfeasible, and + should write in Cut the set of names of nodes in the unfeasible cut (note + that node names depend on USENAME0), Inf<Index>()-terminated, returning the + deficit of the cut (which allows to distinguish which of the two cases + above hold). In general, no special properties can be expected from the + returned cut, but solvers should be able to provide e.g. "small" cuts. + + However, not all solvers may be (easily) capable of providing this + information; thus, returning 0 (no cut) is allowed, as in the base class + implementation, to signify that this information is not available. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Getting unboundedness certificate + @{ */ + + virtual Index MCFGetUnbCycl( Index_Set Pred , Index_Set ArcPred ) + { + return( Inf<Index>() ); + } + +/**< Return an unboundedness certificate. In an unbounded MCF problem, + unboundedness can always be reduced to the existence of a directed cycle + with negative cost and all arcs having infinite capacity. + When detecting unboundedness, MCF solvers are typically capable to provide + one such cycle. This information can be useful, and this method is provided + for getting it. It can be called only if MCFGetStatus() == kUnbounded, and + writes in Pred[] and ArcPred[], respectively, the node and arc predecessor + function of the cycle. That is, if node `i' belongs to the cycle then + `Pred[ i ]' contains the name of the predecessor of `j' of `i' in the cycle + (note that node names depend on USENAME0), and `ArcPred[ i ]' contains the + index of the arc joining the two (note that in general there may be + multiple copies of each arc). Entries of the vectors for nodes not + belonging to the cycle are in principle undefined, and the name of one node + belonging to the cycle is returned by the method. + Note that if there are multiple cycles with negative costs this method + will return just one of them (finding the cycle with most negative cost + is an NO-hard problem), although solvers should be able to produce cycles + with "large negative" cost. + + However, not all solvers may be (easily) capable of providing this + information; thus, returning Inf<Index>() is allowed, as in the base class + implementation, to signify that this information is not available. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Saving/restoring the state of the solver + @{ */ + + virtual MCFStatePtr MCFGetState( void ) + { + return( NULL ); + } + + /**< Save the state of the MCF solver. The MCFClass interface supports the + notion of saving and restoring the state of the MCF solver, such as the + current/optimal basis in a simplex solver. The "empty" class MCFState is + defined as a placeholder for state descriptions. + + MCFGetState() creates and returns a pointer to an object of (a proper + derived class of) class MCFState which describes the current state of the + MCF solver. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void MCFPutState( MCFStatePtr S ) {} + +/**< Restore the solver to the state in which it was when the state `S' was + created with MCFGetState() [see above]. + + The typical use of this method is the following: a MCF problem is solved + and the "optimal state" is set aside. Then the problem changes and it is + re-solved. Then, the problem has to be changed again to a form that is + close to the original one but rather different from the second one (think + of a long backtracking in a Branch & Bound) to which the current "state" + referes. Then, the old optimal state can be expected to provide a better + starting point for reoptimization [see ReOptimize() below]. + + Note, however, that the state is only relative to the optimization process, + i.e., this operation is meaningless if the data of the problem has changed + in the meantime. So, if a state has to be used for speeding up + reoptimization, the following has to be done: + + - first, the data of the solver is brought back to *exactly* the same as + it was at the moment where the state `S' was created (prior than this + operation a ReOptimize( false ) call is probably advisable); + + - then, MCFPutState() is called (and ReOptimize( true ) is called); + + - only afterwards the data of the problem is changed to the final state + and the problem is solved. + + A "put state" operation does not "deplete" the state, which can therefore + be used more than once. Indeed, a state is constructed inside the solver + for each call to MCFGetState(), but the solver never deletes statuses; + this has to be done on the outside when they are no longer needed (the + solver must be "resistent" to deletion of the state at any moment). + + Since not all the MCF solvers reoptimize (efficiently enough to make these + operations worth), an "empty" implementation that does nothing is provided + by the base class. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Time the code + @{ */ + + void TimeMCF( double &t_us , double &t_ss ) + { + t_us = t_ss = 0; + if( MCFt ) + MCFt->Read( t_us , t_ss ); + } + +/**< Time the code. If called within any of the methods of the class that are + "actively timed" (this depends on the subclasses), this method returns the + user and sistem time (in seconds) since the start of that method. If + methods that are actively timed call other methods that are actively + timed, TimeMCF() returns the (...) time since the beginning of the + *outer* actively timed method. If called outside of any actively timed + method, this method returns the (...) time spent in all the previous + executions of all the actively timed methods of the class. + + Implementing the proper calls to MCFt->Start() and MCFt->Stop() is due to + derived classes; these should at least be placed at the beginning and at + the end, respectively, of SolveMCF() and presumably the Chg***() methods, + that is, at least these methods should be "actively timed". */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + double TimeMCF( void ) + { + return( MCFt ? MCFt->Read() : 0 ); + } + +/**< Like TimeMCF(double,double) [see above], but returns the total time. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Check the solutions + @{ */ + + inline void CheckPSol( void ); + +/**< Check that the primal solution returned by the solver is primal feasible. + (to within the tolerances set by SetPar(kEps****) [see above], if any). Also, + check that the objective function value is correct. + + This method is implemented by the base class, using the above methods + for collecting the solutions and the methods of the next section for + reading the data of the problem; as such, they will work for any derived + class that properly implements all these methods. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + inline void CheckDSol( void ); + +/**< Check that the dual solution returned by the solver is dual feasible. + (to within the tolerances set by SetPar(kEps****) [see above], if any). Also, + check that the objective function value is correct. + + This method is implemented by the base class, using the above methods + for collecting the solutions and the methods of the next section for + reading the data of the problem; as such, they will work for any derived + class that properly implements all these methods. */ + +/*@} -----------------------------------------------------------------------*/ +/*-------------- METHODS FOR READING THE DATA OF THE PROBLEM ---------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Reading graph size + @{ */ + + inline Index MCFnmax( void ) + { + return( nmax ); + } + +/**< Return the maximum number of nodes for this instance of MCFClass. + The implementation of the method in the base class returns the protected + fields \c nmax, which is provided for derived classes to hold this + information. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + inline Index MCFmmax( void ) + { + return( mmax ); + } + +/**< Return the maximum number of arcs for this instance of MCFClass. + The implementation of the method in the base class returns the protected + fields \c mmax, which is provided for derived classes to hold this + information. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + inline Index MCFn( void ) + { + return( n ); + } + +/**< Return the number of nodes in the current graph. + The implementation of the method in the base class returns the protected + fields \c n, which is provided for derived classes to hold this + information. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + inline Index MCFm( void ) + { + return( m ); + } + +/**< Return the number of arcs in the current graph. + The implementation of the method in the base class returns the protected + fields \c m, which is provided for derived classes to hold this + information. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading graph topology + @{ */ + + virtual void MCFArcs( Index_Set Startv , Index_Set Endv , + cIndex_Set nms = NULL , cIndex strt = 0 , + Index stp = Inf<Index>() ) = 0; + +/**< Write the starting (tail) and ending (head) nodes of the arcs in Startv[] + and Endv[]. If nms == NULL, then the information relative to all arcs is + written into Startv[] and Endv[], otherwise Startv[ i ] and Endv[ i ] + contain the information relative to arc nms[ i ] (nms[] must be + Inf<Index>()-terminated). + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to arcs which are *both* in nms and whose index is in + the correct range are returned. + + Startv or Endv can be NULL, meaning that only the other information is + required. + + @note If USENAME0 == 0 then the returned node names will be in the range + 1 .. n, while if USENAME0 == 1 the returned node names will be in + the range 0 .. n - 1. + + @note If the graph is "dynamic", be careful to use MCFn() e MCFm() to + properly choose the dimension of nodes and arcs arrays. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual Index MCFSNde( cIndex i ) = 0; + +/**< Return the starting (tail) node of the arc `i'. + + @note If USENAME0 == 0 then the returned node names will be in the range + 1 .. n, while if USENAME0 == 1 the returned node names will be in + the range 0 .. n - 1. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual Index MCFENde( cIndex i ) = 0; + +/**< Return the ending (head) node of the arc `i'. + + @note If USENAME0 == 0 then the returned node names will be in the range + 1 .. n, while if USENAME0 == 1 the returned node names will be in + the range 0 .. n - 1. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cIndex_Set MCFSNdes( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal vector containing the starting + (tail) nodes for each arc. Since this may *not always be available*, + depending on the implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cIndex_Set MCFENdes( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal vector containing the ending + (head) nodes for each arc. Since this may *not always be available*, + depending on the implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading arc costs + @{ */ + + virtual void MCFCosts( CRow Costv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Write the arc costs into Costv[]. If nms == NULL, then all the costs are + written, otherwise Costv[ i ] contains the information relative to arc + nms[ i ] (nms[] must be Inf<Index>()-terminated). + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to arcs which are *both* in nms and whose index is in + the correct range are returned. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual CNumber MCFCost( cIndex i ) = 0; + +/**< Return the cost of the i-th arc. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cCRow MCFCosts( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal vector containing the arc + costs. Since this may *not always be available*, depending on the + implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void MCFQCoef( CRow Qv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) + +/**< Write the quadratic coefficients of the arc costs into Qv[]. If + nms == NULL, then all the costs are written, otherwise Costv[ i ] contains + the information relative to arc nms[ i ] (nms[] must be + Inf<Index>()-terminated). + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to arcs which are *both* in nms and whose index is in + the correct range are returned. + + Note that the method is *not* pure virtual: an implementation is provided + for "pure linear" MCF solvers that only work with all zero quadratic + coefficients. */ + { + if( nms ) { + while( *nms < strt ) + nms++; + + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Qv++) = 0; + } + else { + if( stp > m ) + stp = m; + + while( stp-- > strt ) + *(Qv++) = 0; + } + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual CNumber MCFQCoef( cIndex i ) + +/**< Return the quadratic coefficients of the cost of the i-th arc. Note that + the method is *not* pure virtual: an implementation is provided for "pure + linear" MCF solvers that only work with all zero quadratic coefficients. */ + { + return( 0 ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cCRow MCFQCoef( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal vector containing the arc + costs. Since this may *not always be available*, depending on the + implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready (such as "pure linear" MCF solvers that only + work with all zero quadratic coefficients) does not need to implement the + method. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading arc capacities + @{ */ + + virtual void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Write the arc capacities into UCapv[]. If nms == NULL, then all the + capacities are written, otherwise UCapv[ i ] contains the information + relative to arc nms[ i ] (nms[] must be Inf<Index>()-terminated). + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to arcs which are *both* in nms and whose index is in + the correct range are returned. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual FNumber MCFUCap( cIndex i ) = 0; + +/**< Return the capacity of the i-th arc. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cFRow MCFUCaps( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal vector containing the arc + capacities. Since this may *not always be available*, depending on the + implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Reading node deficits + @{ */ + + virtual void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Write the node deficits into Dfctv[]. If nms == NULL, then all the + defcits are written, otherwise Dfctvv[ i ] contains the information + relative to node nms[ i ] (nms[] must be Inf<Index>()-terminated). + + The parameters `strt' and `stp' allow to restrict the output of the method + to all and only the nodes `i' with strt <= i < min( MCFn() , stp ). `strt' + and `stp' work in "&&" with nms; that is, if nms != NULL then only the + values corresponding to nodes which are *both* in nms and whose index is in + the correct range are returned. + + @note Here node "names" (strt and stp, those contained in nms[] or `i' + in MCFDfct()) go from 0 to n - 1, regardless to the value of + USENAME0; hence, if USENAME0 == 0 then the first node is "named 1", + but its deficit is returned by MCFDfcts( Dfctv , NULL , 0 , 1 ). */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual FNumber MCFDfct( cIndex i ) = 0; + +/**< Return the deficit of the i-th node. + + @note Here node "names" go from 0 to n - 1, regardless to the value of + USENAME0; hence, if USENAME0 == 0 then the first node is "named 1", + but its deficit is returned by MCFDfct( 0 ). */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual cFRow MCFDfcts( void ) + { + return( NULL ); + } + +/**< Return a read-only pointer to an internal vector containing the node + deficits. Since this may *not always be available*, depending on the + implementation, this method can (uniformly) return NULL. + This is done by the base class already, so a derived class that does not + have the information ready does not need to implement the method. + + @note Here node "names" go from 0 to n - 1, regardless to the value of + USENAME0; hence, if USENAME0 == 0 then the first node is "named 1", + but its deficit is contained in MCFDfcts()[ 0 ]. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Write problem to file + @{ */ + + virtual inline void WriteMCF( ostream &oStrm , int frmt = 0 ); + +/**< Write the current MCF problem to an ostream. This may be useful e.g. for + debugging purposes. + + The base MCFClass class provides output in two different formats, depending + on the value of the parameter frmt: + + - kDimacs the problem is written in DIMACS standard format, read by most + MCF codes available; + + - kMPS the problem is written in the "modern version" (tab-separated) + of the MPS format, read by most LP/MIP solvers; + + - kFWMPS the problem is written in the "old version" (fixed width + fields) of the MPS format; this is read by most LP/MIP + solvers, but some codes still require the old format. + + The implementation of WriteMCF() in the base class uses all the above + methods for reading the data; as such it will work for any derived class + that properly implements this part of the interface, but it may not be + very efficient. Thus, the method is virtual to allow the derived classes + to either re-implement WriteMCF() for the above two formats in a more + efficient way, and/or to extend it to support other solver-specific + formats. + + @note None of the above two formats supports quadratic MCFs, so if + nonzero quadratic coefficients are present, they are just ignored. + */ + +/*@} -----------------------------------------------------------------------*/ +/*------------- METHODS FOR ADDING / REMOVING / CHANGING DATA --------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Changing the costs + @{ */ + + virtual void ChgCosts( cCRow NCost , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Change the arc costs. In particular, change the costs that are: + + - listed in into the vector of indices `nms' (ordered in increasing + sense and Inf<Index>()-terminated), + + - *and* whose name belongs to the interval [`strt', `stp'). + + That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th cost will be + changed to NCost[ i ]. If nms == NULL (as the default), *all* the + entries in the given range will be changed; if stp > MCFm(), then the + smaller bound is used. + + @note changing the costs of arcs that *do not exist* is *not allowed*; + only arcs which have not been closed/deleted [see CloseArc() / + DelArc() below and LoadNet() above about C_INF costs] can be + touched with these methods. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void ChgCost( Index arc , cCNumber NCost ) = 0; + +/**< Change the cost of the i-th arc. + + @note changing the costs of arcs that *do not exist* is *not allowed*; + only arcs which have not been closed/deleted [see CloseArc() / + DelArc() below and LoadNet() above about C_INF costs] can be + touched with these methods. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) + +/**< Change the quadratic coefficients of the arc costs. In particular, + change the coefficients that are: + + - listed in into the vector of indices `nms' (ordered in increasing + sense and Inf<Index>()-terminated), + + - *and* whose name belongs to the interval [`strt', `stp'). + + That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th cost will be + changed to NCost[ i ]. If nms == NULL (as the default), *all* the + entries in the given range will be changed; if stp > MCFm(), then the + smaller bound is used. If NQCoef == NULL, all the specified coefficients + are set to zero. + + Note that the method is *not* pure virtual: an implementation is provided + for "pure linear" MCF solvers that only work with all zero quadratic + coefficients. + + @note changing the costs of arcs that *do not exist* is *not allowed*; + only arcs which have not been closed/deleted [see CloseArc() / + DelArc() below and LoadNet() above about C_INF costs] can be + touched with these methods. */ + { + assert(NQCoef); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void ChgQCoef( Index arc , cCNumber NQCoef ) + +/**< Change the quadratic coefficient of the cost of the i-th arc. + + Note that the method is *not* pure virtual: an implementation is provided + for "pure linear" MCF solvers that only work with all zero quadratic + coefficients. + + @note changing the costs of arcs that *do not exist* is *not allowed*; + only arcs which have not been closed/deleted [see CloseArc() / + DelArc() below and LoadNet() above about C_INF costs] can be + touched with these methods. */ + { + assert(NQCoef); + } + +/*@} -----------------------------------------------------------------------*/ +/** @name Changing the capacities + @{ */ + + virtual void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Change the arc capacities. In particular, change the capacities that are: + + - listed in into the vector of indices `nms' (ordered in increasing sense + and Inf<Index>()-terminated), + + - *and* whose name belongs to the interval [`strt', `stp'). + + That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th capacity will be + changed to NCap[ i ]. If nms == NULL (as the default), *all* the + entries in the given range will be changed; if stp > MCFm(), then the + smaller bound is used. + + @note changing the capacities of arcs that *do not exist* is *not allowed*; + only arcs that have not been closed/deleted [see CloseArc() / + DelArc() below and LoadNet() above about C_INF costs] can be + touched with these methods. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void ChgUCap( Index arc , cFNumber NCap ) = 0; + +/**< Change the capacity of the i-th arc. + + @note changing the capacities of arcs that *do not exist* is *not allowed*; + only arcs that have not been closed/deleted [see CloseArc() / + DelArc() below and LoadNet() above about C_INF costs] can be + touched with these methods. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Changing the deficits + @{ */ + + virtual void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) = 0; + +/**< Change the node deficits. In particular, change the deficits that are: + + - listed in into the vector of indices `nms' (ordered in increasing sense + and Inf<Index>()-terminated), + + - *and* whose name belongs to the interval [`strt', `stp'). + + That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th deficit will be + changed to NDfct[ i ]. If nms == NULL (as the default), *all* the + entries in the given range will be changed; if stp > MCFn(), then the + smaller bound is used. + + Note that, in ChgDfcts(), node "names" (strt, stp or those contained + in nms[]) go from 0 to n - 1, regardless to the value of USENAME0; hence, + if USENAME0 == 0 then the first node is "named 1", but its deficit can be + changed e.g. with ChgDfcts( &new_deficit , NULL , 0 , 1 ). + + @note changing the capacities of nodes that *do not exist* is *not + allowed*; only nodes that have not been deleted [see DelNode() + below] can be touched with these methods. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void ChgDfct( Index node , cFNumber NDfct ) = 0; + +/**< Change the deficit of the i-th node. + + Note that, in ChgDfct[s](), node "names" (i, strt/ stp or those contained + in nms[]) go from 0 to n - 1, regardless to the value of USENAME0; hence, + if USENAME0 == 0 then the first node is "named 1", but its deficit can be + changed e.g. with ChgDfcts( 0 , new_deficit ). + + @note changing the capacities of nodes that *do not exist* is *not + allowed*; only nodes that have not been deleted [see DelNode() + below] can be touched with these methods. */ + +/*@} -----------------------------------------------------------------------*/ +/** @name Changing graph topology + @{ */ + + virtual void CloseArc( cIndex name ) = 0; + +/**< "Close" the arc `name'. Although all the associated information (name, + cost, capacity, end and start node) is kept, the arc is removed from the + problem until OpenArc( i ) [see below] is called. + + "closed" arcs always have 0 flow, but are otherwise counted as any other + arc; for instance, MCFm() does *not* decrease as an effect of a call to + CloseArc(). How this closure is implemented is solver-specific. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual bool IsClosedArc( cIndex name ) = 0; + +/**< IsClosedArc() returns true if and only if the arc `name' is closed. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void DelNode( cIndex name ) = 0; + +/**< Delete the node `name'. + + For any value of `name', all incident arcs to that node are closed [see + CloseArc() above] (*not* Deleted, see DelArc() below) and the deficit is + set to zero. + + Il furthermore `name' is the last node, the number of nodes as reported by + MCFn() is reduced by at least one, until the n-th node is not a deleted + one. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void OpenArc( cIndex name ) = 0; + +/**< Restore the previously closed arc `name'. It is an error to open an arc + that has not been previously closed. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual Index AddNode( cFNumber aDfct ) = 0; + +/**< Add a new node with deficit aDfct, returning its name. Inf<Index>() is + returned if there is no room for a new node. Remember that the node names + are either { 0 .. nmax - 1 } or { 1 .. nmax }, depending on the value of + USENAME0. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void ChangeArc( cIndex name , cIndex nSN = Inf<Index>() , + cIndex nEN = Inf<Index>() ) = 0; + +/**< Change the starting and/or ending node of arc `name' to nSN and nEN. + Each parameter being Inf<Index>() means to leave the previous starting or + ending node untouched. When this method is called `name' can be either the + name of a "normal" arc or that of a "closed" arc [see CloseArc() above]: + in the latter case, at the end of ChangeArc() the arc is *still closed*, + and it remains so until OpenArc( name ) [see above] is called. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual void DelArc( cIndex name ) = 0; + +/**< Delete the arc `name'. Unlike "closed" arcs, all the information + associated with a deleted arc is lost and `name' is made available as a + name for new arcs to be created with AddArc() [see below]. + + Il furthermore `name' is the last arc, the number of arcs as reported by + MCFm() is reduced by at least one, until the m-th arc is not a deleted + one. Otherwise, the flow on the arc is always ensured to be 0. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual bool IsDeletedArc( cIndex name ) = 0; + +/**< Return true if and only if the arc `name' is deleted. It should only + be called with name < MCFm(), as every other arc is deleted by + definition. */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + virtual Index AddArc( cIndex Start , cIndex End , cFNumber aU , + cCNumber aC ) = 0; + +/**< Add the new arc ( Start , End ) with cost aC and capacity aU, + returning its name. Inf<Index>() is returned if there is no room for a new + arc. Remember that arc names go from 0 to mmax - 1. */ + +/*@} -----------------------------------------------------------------------*/ +/*------------------------------ DESTRUCTOR --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Destructor + @{ */ + + virtual ~MCFClass() + { + delete MCFt; + } + +/**< Destructor of the class. The implementation in the base class only + deletes the MCFt field. It is virtual, as it should be. */ + +/*@} -----------------------------------------------------------------------*/ +/*-------------------- PROTECTED PART OF THE CLASS -------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-- --*/ +/*-- The following fields are believed to be general enough to make it --*/ +/*-- worth adding them to the abstract base class. --*/ +/*-- --*/ +/*--------------------------------------------------------------------------*/ + + protected: + +/*--------------------------------------------------------------------------*/ +/*--------------------------- PROTECTED METHODS ----------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-------------------------- MANAGING COMPARISONS --------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @name Managing comparisons. + The following methods are provided for making it easier to perform + comparisons, with and without tolerances. + @{ */ + +/*--------------------------------------------------------------------------*/ +/** true if flow x is equal to zero (possibly considering tolerances). */ + + template<class T> + inline bool ETZ( T x , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x == 0 ); + else + return( (x <= eps ) && ( x >= -eps ) ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is greater than zero (possibly considering tolerances). */ + + template<class T> + inline bool GTZ( T x , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x > 0 ); + else + return( x > eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is greater than or equal to zero (possibly considering + tolerances). */ + + template<class T> + inline bool GEZ( T x , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x >= 0 ); + else + return( x >= - eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is less than zero (possibly considering tolerances). */ + + template<class T> + inline bool LTZ( T x , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x < 0 ); + else + return( x < - eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is less than or equal to zero (possibly considering + tolerances). */ + + template<class T> + inline bool LEZ( T x , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x <= 0 ); + else + return( x <= eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is greater than flow y (possibly considering tolerances). + */ + + template<class T> + inline bool GT( T x , T y , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x > y ); + else + return( x > y + eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is less than flow y (possibly considering tolerances). */ + + template<class T> + inline bool LT( T x , T y , const T eps ) + { + if( numeric_limits<T>::is_integer ) + return( x < y ); + else + return( x < y - eps ); + } + +/*@} -----------------------------------------------------------------------*/ +/*---------------------- PROTECTED DATA STRUCTURES -------------------------*/ +/*--------------------------------------------------------------------------*/ + + Index n; ///< total number of nodes + Index nmax; ///< maximum number of nodes + + Index m; ///< total number of arcs + Index mmax; ///< maximum number of arcs + + int status; ///< return status, see the comments to MCFGetStatus() above. + ///< Note that the variable is defined int to allow derived + ///< classes to return their own specialized status codes + bool Senstv; ///< true <=> the latest optimal solution should be exploited + + OPTtimers *MCFt; ///< timer for performances evaluation + + FNumber EpsFlw; ///< precision for comparing arc flows / capacities + FNumber EpsDfct; ///< precision for comparing node deficits + CNumber EpsCst; ///< precision for comparing arc costs + + double MaxTime; ///< max time (in seconds) in which MCF Solver can find + ///< an optimal solution (0 = no limits) + int MaxIter; ///< max number of iterations in which MCF Solver can find + ///< an optimal solution (0 = no limits) + +/*--------------------------------------------------------------------------*/ + + }; // end( class MCFClass ) + +/*@} -----------------------------------------------------------------------*/ +/*------------------- inline methods implementation ------------------------*/ +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::LoadDMX( istream &DMXs , bool IsQuad ) +{ + // read first non-comment line - - - - - - - - - - - - - - - - - - - - - - - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + char c; + for(;;) { + if( ! ( DMXs >> c ) ) + throw( MCFException( "LoadDMX: error reading the input stream" ) ); + + if( c != 'c' ) // if it's not a comment + break; + + do // skip the entire line + if( ! DMXs.get( c ) ) + throw( MCFException( "LoadDMX: error reading the input stream" ) ); + while( c != '\n' ); + } + + if( c != 'p' ) + throw( MCFException( "LoadDMX: format error in the input stream" ) ); + + char buf[ 3 ]; // free space + DMXs >> buf; // skip (in has to be "min") + + Index tn; + if( ! ( DMXs >> tn ) ) + throw( MCFException( "LoadDMX: error reading number of nodes" ) ); + + Index tm; + if( ! ( DMXs >> tm ) ) + throw( MCFException( "LoadDMX: error reading number of arcs" ) ); + + // allocate memory - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + FRow tU = new FNumber[ tm ]; // arc upper capacities + FRow tDfct = new FNumber[ tn ]; // node deficits + Index_Set tStartn = new Index[ tm ]; // arc start nodes + Index_Set tEndn = new Index[ tm ]; // arc end nodes + CRow tC = new CNumber[ tm ]; // arc costs + CRow tQ = NULL; + if( IsQuad ) + tQ = new CNumber[ tm ]; // arc quadratic costs + + + for( Index i = 0 ; i < tn ; ) // all deficits are 0 + tDfct[ i++ ] = 0; // unless otherwise stated + + // read problem data - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + Index i = 0; // arc counter + for(;;) { + if( ! ( DMXs >> c ) ) + break; + + switch( c ) { + case( 'c' ): // comment line- - - - - - - - - - - - - - - - - - - - - - - + do + if( ! DMXs.get( c ) ) + throw( MCFException( "LoadDMX: error reading the input stream" ) ); + while( c != '\n' ); + break; + + case( 'n' ): // description of a node - - - - - - - - - - - - - - - - - - + Index j; + if( ! ( DMXs >> j ) ) + throw( MCFException( "LoadDMX: error reading node name" ) ); + + if( ( j < 1 ) || ( j > tn ) ) + throw( MCFException( "LoadDMX: invalid node name" ) ); + + FNumber Dfctj; + if( ! ( DMXs >> Dfctj ) ) + throw( MCFException( "LoadDMX: error reading deficit" ) ); + + tDfct[ j - 1 ] -= Dfctj; + break; + + case( 'a' ): // description of an arc - - - - - - - - - - - - - - - - - - + if( i == tm ) + throw( MCFException( "LoadDMX: too many arc descriptors" ) ); + + if( ! ( DMXs >> tStartn[ i ] ) ) + throw( MCFException( "LoadDMX: error reading start node" ) ); + + if( ( tStartn[ i ] < 1 ) || ( tStartn[ i ] > tn ) ) + throw( MCFException( "LoadDMX: invalid start node" ) ); + + if( ! ( DMXs >> tEndn[ i ] ) ) + throw( MCFException( "LoadDMX: error reading end node" ) ); + + if( ( tEndn[ i ] < 1 ) || ( tEndn[ i ] > tn ) ) + throw( MCFException( "LoadDMX: invalid end node" ) ); + + if( tStartn[ i ] == tEndn[ i ] ) + throw( MCFException( "LoadDMX: self-loops not permitted" ) ); + + FNumber LB; + if( ! ( DMXs >> LB ) ) + throw( MCFException( "LoadDMX: error reading lower bound" ) ); + + if( ! ( DMXs >> tU[ i ] ) ) + throw( MCFException( "LoadDMX: error reading upper bound" ) ); + + if( ! ( DMXs >> tC[ i ] ) ) + throw( MCFException( "LoadDMX: error reading arc cost" ) ); + + if( tQ ) { + if( ! ( DMXs >> tQ[ i ] ) ) + throw( MCFException( "LoadDMX: error reading arc quadratic cost" ) ); + + if( tQ[ i ] < 0 ) + throw( MCFException( "LoadDMX: negative arc quadratic cost" ) ); + } + + if( tU[ i ] < LB ) + throw( MCFException( "LoadDMX: lower bound > upper bound" ) ); + + tU[ i ] -= LB; + tDfct[ tStartn[ i ] - 1 ] += LB; + tDfct[ tEndn[ i ] - 1 ] -= LB; + #if( USENAME0 ) + tStartn[ i ]--; // in the DIMACS format, node names start from 1 + tEndn[ i ]--; + #endif + i++; + break; + + default: // invalid code- - - - - - - - - - - - - - - - - - - - - - - - - + throw( MCFException( "LoadDMX: invalid code" ) ); + + } // end( switch( c ) ) + } // end( for( ever ) ) + + if( i < tm ) + throw( MCFException( "LoadDMX: too few arc descriptors" ) ); + + // call LoadNet- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + LoadNet( tn , tm , tn , tm , tU , tC , tDfct , tStartn , tEndn ); + + // then pass quadratic costs, if any + + if( tQ ) + ChgQCoef( tQ ); + + // delete the original data structures - - - - - - - - - - - - - - - - - - - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + delete[] tQ; + delete[] tC; + delete[] tEndn; + delete[] tStartn; + delete[] tDfct; + delete[] tU; + + } // end( MCFClass::LoadDMX ) + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::SetPar( int par , int val ) +{ + switch( par ) { + case( kMaxIter ): + MaxIter = val; + break; + + case( kReopt ): + Senstv = (val == kYes); + break; + + default: + throw( MCFException( "Error using SetPar: unknown parameter" ) ); + } + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::SetPar( int par, double val ) +{ + switch( par ) { + case( kEpsFlw ): + if( EpsFlw != FNumber( val ) ) { + EpsFlw = FNumber( val ); + EpsDfct = EpsFlw * ( nmax ? nmax : 100 ); + status = kUnSolved; + } + break; + + case( kEpsDfct ): + if( EpsDfct != FNumber( val ) ) { + EpsDfct = FNumber( val ); + status = kUnSolved; + } + break; + + case( kEpsCst ): + if( EpsCst != CNumber( val ) ) { + EpsCst = CNumber( val ); + status = kUnSolved; + } + break; + + case( kMaxTime ): + MaxTime = val; + break; + + default: + throw( MCFException( "Error using SetPar: unknown parameter" ) ); + } + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::GetPar( int par, int &val ) +{ + switch( par ) { + case( kMaxIter ): + val = MaxIter; + break; + + case( kReopt ): + if( Senstv ) val = kYes; + else val = kNo; + break; + + default: + throw( MCFException( "Error using GetPar: unknown parameter" ) ); + } + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::GetPar( int par , double &val ) +{ + switch( par ) { + case( kEpsFlw ): + val = double( EpsFlw ); + break; + + case( kEpsDfct ): + val = double( EpsDfct ); + break; + + case( kEpsCst ): + val = double( EpsCst ); + break; + + case( kMaxTime ): + val = MaxTime; + break; + + default: + throw( MCFException( "Error using GetPar: unknown parameter" ) ); + } + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::CheckPSol( void ) +{ + FRow tB = new FNumber[ MCFn() ]; + MCFDfcts( tB ); + #if( ! USENAME0 ) + tB--; + #endif + + FRow tX = new FNumber[ MCFm() ]; + MCFGetX( tX ); + FONumber CX = 0; + for( Index i = 0 ; i < MCFm() ; i++ ) + if( GTZ( tX[ i ] , EpsFlw ) ) { + if( IsClosedArc( i ) ) + throw( MCFException( "Closed arc with nonzero flow (CheckPSol)" ) ); + + if( GT( tX[ i ] , MCFUCap( i ) , EpsFlw ) ) + throw( MCFException( "Arc flow exceeds capacity (CheckPSol)" ) ); + + if( LTZ( tX[ i ] , EpsFlw ) ) + throw( MCFException( "Arc flow is negative (CheckPSol)" ) ); + + CX += ( FONumber( MCFCost( i ) ) + + FONumber( MCFQCoef( i ) ) * FONumber( tX[ i ] ) / 2 ) + * FONumber( tX[ i ] ); + tB[ MCFSNde( i ) ] += tX[ i ]; + tB[ MCFENde( i ) ] -= tX[ i ]; + } + delete[] tX; + #if( ! USENAME0 ) + tB++; + #endif + + for( Index i = 0 ; i < MCFn() ; i++ ) + if( ! ETZ( tB[ i ] , EpsDfct ) ) + throw( MCFException( "Node is not balanced (CheckPSol)" ) ); + + delete[] tB; + CX -= MCFGetFO(); + if( ( CX >= 0 ? CX : - CX ) > EpsCst * MCFn() ) + throw( MCFException( "Objective function value is wrong (CheckPSol)" ) ); + + } // end( CheckPSol ) + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::CheckDSol( void ) +{ + CRow tPi = new CNumber[ MCFn() ]; + MCFGetPi( tPi ); + + FONumber BY = 0; + for( Index i = 0 ; i < MCFn() ; i++ ) + BY += FONumber( tPi[ i ] ) * FONumber( MCFDfct( i ) ); + + CRow tRC = new CNumber[ MCFm() ]; + MCFGetRC( tRC ); + + FRow tX = new FNumber[ MCFm() ]; + MCFGetX( tX ); + + #if( ! USENAME0 ) + tPi--; + #endif + + FONumber QdrtcCmpnt = 0; // the quadratic part of the objective + for( Index i = 0 ; i < MCFm() ; i++ ) { + if( IsClosedArc( i ) ) + continue; + + cFONumber tXi = FONumber( tX[ i ] ); + cCNumber QCi = CNumber( MCFQCoef( i ) * tXi ); + QdrtcCmpnt += QCi * tXi; + CNumber RCi = MCFCost( i ) + QCi + + tPi[ MCFSNde( i ) ] - tPi[ MCFENde( i ) ]; + RCi -= tRC[ i ]; + + if( ! ETZ( RCi , EpsCst ) ) + throw( MCFException( "Reduced Cost value is wrong (CheckDSol)" ) ); + + if( LTZ( tRC[ i ] , EpsCst ) ) { + BY += FONumber( tRC[ i ] ) * FONumber( MCFUCap( i ) ); + + if( LT( tX[ i ] , MCFUCap( i ) , EpsFlw ) ) + throw( MCFException( "Csomplementary Slackness violated (CheckDSol)" ) ); + } + else + if( GTZ( tRC[ i ] , EpsCst ) && GTZ( tX[ i ] , EpsFlw ) ) + throw( MCFException( "Complementary Slackness violated (CheckDSol)" ) ); + + } // end( for( i ) ) + + delete[] tX; + delete[] tRC; + #if( ! USENAME0 ) + tPi++; + #endif + delete[] tPi; + + BY -= ( MCFGetFO() + QdrtcCmpnt / 2 ); + if( ( BY >= 0 ? BY : - BY ) > EpsCst * MCFn() ) + throw( MCFException( "Objective function value is wrong (CheckDSol)" ) ); + + } // end( CheckDSol ) + +/*--------------------------------------------------------------------------*/ + +inline void MCFClass::WriteMCF( ostream &oStrm , int frmt ) +{ + if( ( ! numeric_limits<FNumber>::is_integer ) || + ( ! numeric_limits<CNumber>::is_integer ) ) + oStrm.precision( 12 ); + + switch( frmt ) { + case( kDimacs ): // DIMACS standard format - - - - - - - - - - - - - - - - + //- - - - - - - - - - - - - - - - - - - - - - - - - - - - + case( kQDimacs ): // ... linear or quadratic + + // compute and output preamble and size- - - - - - - - - - - - - - - - - - + + oStrm << "p min " << MCFn() << " "; + { + Index tm = MCFm(); + for( Index i = MCFm() ; i-- ; ) + if( IsClosedArc( i ) || IsDeletedArc( i ) ) + tm--; + + oStrm << tm << endl; + } + + // output arc information- - - - - - - - - - - - - - - - - - - - - - - - - + + for( Index i = 0 ; i < MCFm() ; i++ ) + if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) { + oStrm << "a\t"; + #if( USENAME0 ) + oStrm << MCFSNde( i ) + 1 << "\t" << MCFENde( i ) + 1 << "\t"; + #else + oStrm << MCFSNde( i ) << "\t" << MCFENde( i ) << "\t"; + #endif + oStrm << "0\t" << MCFUCap( i ) << "\t" << MCFCost( i ); + + if( frmt == kQDimacs ) + oStrm << "\t" << MCFQCoef( i ); + + oStrm << endl; + } + + // output node information - - - - - - - - - - - - - - - - - - - - - - - - + + for( Index i = 0 ; i < MCFn() ; ) { + cFNumber Dfcti = MCFDfct( i++ ); + if( Dfcti ) + oStrm << "n\t" << i << "\t" << - Dfcti << endl; + } + + break; + + case( kMPS ): // MPS format(s)- - - - - - - - - - - - - - - - - - - - - + case( kFWMPS ): //- - - - - - - - - - - - - - - - - - - - - - - - - - - - + + // writing problem name- - - - - - - - - - - - - - - - - - - - - - - - - - + + oStrm << "NAME MCF" << endl; + + // writing ROWS section: flow balance constraints- - - - - - - - - - - - - + + oStrm << "ROWS" << endl << " N obj" << endl; + + for( Index i = 0 ; i < MCFn() ; ) + oStrm << " E c" << i++ << endl; + + // writing COLUMNS section - - - - - - - - - - - - - - - - - - - - - - - - + + oStrm << "COLUMNS" << endl; + + if( frmt == kMPS ) { // "modern" MPS + for( Index i = 0 ; i < MCFm() ; i++ ) + if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) { + oStrm << " x" << i << "\tobj\t" << MCFCost( i ) << "\tc"; + #if( USENAME0 ) + oStrm << MCFSNde( i ) << "\t-1" << endl << " x" << i << "\tc" + << MCFENde( i ) << "\t1" << endl; + #else + oStrm << MCFSNde( i ) - 1 << "\t-1" << endl << " x" << i << "\tc" + << MCFENde( i ) - 1 << "\t1" << endl; + #endif + } + } + else // "old" MPS + for( Index i = 0 ; i < MCFm() ; i++ ) { + ostringstream myField; + myField << "x" << i << ends; + oStrm << " " << setw( 8 ) << left << myField.str() + << " " << setw( 8 ) << left << "obj" + << " " << setw( 12 ) << left << MCFCost( i ); + + myField.seekp( 0 ); + myField << "c" << MCFSNde( i ) - ( 1 - USENAME0 ) << ends; + + oStrm << " " << setw( 8 ) << left << myField.str() + << " " << setw( 12 ) << left << -1 << endl; + + myField.seekp( 0 ); + myField << "x" << i << ends; + + oStrm << " " << setw( 8 ) << left << myField.str(); + + myField.seekp( 0 ); + myField << "c" << MCFENde( i ) - ( 1 - USENAME0 ) << ends; + + oStrm << " " << setw( 8 ) << left << myField.str() << endl; + } + + // writing RHS section: flow balance constraints - - - - - - - - - - - - - + + oStrm << "RHS" << endl; + + if( frmt == kMPS ) // "modern" MPS + for( Index i = 0 ; i < MCFn() ; i++ ) + oStrm << " rhs\tc" << i << "\t" << MCFDfct( i ) << endl; + else // "old" MPS + for( Index i = 0 ; i < MCFn() ; i++ ) { + ostringstream myField; + oStrm << " " << setw( 8 ) << left << "rhs"; + myField << "c" << i << ends; + + oStrm << " " << setw( 8 ) << left << myField.str() + << " " << setw( 12 ) << left << MCFDfct( i ) << endl; + } + + // writing BOUNDS section- - - - - - - - - - - - - - - - - - - - - - - - - + + oStrm << "BOUNDS" << endl; + + if( frmt == kMPS ) { // "modern" MPS + for( Index i = 0 ; i < MCFm() ; i++ ) + if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) + oStrm << " UP BOUND\tx" << i << "\t" << MCFUCap( i ) << endl; + } + else // "old" MPS + for( Index i = 0 ; i < MCFm() ; i++ ) + if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) { + ostringstream myField; + oStrm << " " << setw( 2 ) << left << "UP" + << " " << setw( 8 ) << left << "BOUND"; + + myField << "x" << i << ends; + + oStrm << " " << setw( 8 ) << left << myField.str() + << " " << setw( 12 ) << left << MCFUCap( i ) << endl; + } + + oStrm << "ENDATA" << endl; + break; + + default: // unknown format - - - - - - - - - - - - - - - - - - - - + //- - - - - - - - - - - - - - - - - - - - - - - - - - - - + oStrm << "Error: unknown format " << frmt << endl; + } + } // end( MCFClass::WriteMCF ) + +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) + }; // end( namespace MCFClass_di_unipi_it ) +#endif + +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#undef throw + +#endif /* MCFClass.h included */ + +/*--------------------------------------------------------------------------*/ +/*------------------------- End File MCFClass.h ----------------------------*/ +/*--------------------------------------------------------------------------*/ diff --git a/impl/MCFSimplex.cpp b/impl/MCFSimplex.cpp new file mode 100644 index 0000000..60c70d1 --- /dev/null +++ b/impl/MCFSimplex.cpp @@ -0,0 +1,4198 @@ +/*--------------------------------------------------------------------------*/ +/*---------------------------- File MCFSimplex.C ---------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-- --*/ +/*-- Linear and Quadratic Min Cost Flow problems solver based on the --*/ +/*-- (primal and dual) simplex algorithm. Conforms to the standard MCF --*/ +/*-- interface defined in MCFClass.h. --*/ +/*-- --*/ +/*-- VERSION 1.00 --*/ +/*-- 29 - 08 - 2011 --*/ +/*-- --*/ +/*-- Implementation: --*/ +/*-- --*/ +/*-- Alessandro Bertolini --*/ +/*-- Antonio Frangioni --*/ +/*-- --*/ +/*-- Operations Research Group --*/ +/*-- Dipartimento di Informatica --*/ +/*-- Universita' di Pisa --*/ +/*-- --*/ +/*-- Copyright (C) 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni --*/ +/*-- --*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------- IMPLEMENTATION -------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +/*--------------------------------------------------------------------------*/ +/*------------------------------ INCLUDES ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#include "MCFSimplex.h" +#include <algorithm> +#include <iostream> + +#include <cstdlib> +#include <ctime> + +/*--------------------------------------------------------------------------*/ +/*-------------------------------- USING -----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) +using namespace MCFClass_di_unipi_it; +#endif + +/*--------------------------------------------------------------------------*/ +/*-------------------------------- MACROS ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#define LIMITATEPRECISION 1 + +/* If LIMITATEPRECISION is 1, in the quadratic case the Primal Simplex accepts + entering arc in base only if the decrease of the o.f. value is bigger than + a fixed thresold (EpsOpt * oldFOValue / n). Otherwise, any strict decrease + in the o.f. value is accepted. */ + +#define UNIPI_PRIMAL_INITIAL_SHOW 0 + +/* If UNIPI_PRIMAL_INITIAL_SHOW = 1, Primal Simplex shows the initial condition + (arcs and nodes) of the network. */ + +#define UNIPI_PRIMAL_ITER_SHOW 0 + +/* If UNIPI_PRIMAL_FINAL_SHOW = x with x > 0, Primal Simplex shows the condition + (arcs and nodes) of the network every x iterations. */ + +#define UNIPI_PRIMAL_FINAL_SHOW 0 + +/* If UNIPI_PRIMAL_FINAL_SHOW = 1, Primal Simplex shows the final condition + (arcs and nodes) of the network. */ + +#define UNIPI_DUAL_INITIAL_SHOW 0 + +/* If UNIPI_DUAL_INITIAL_SHOW = 1, Dual Simplex shows the initial condition + (arcs and nodes) of the network. */ + +#define UNIPI_DUAL_ITER_SHOW 0 + +/* If UNIPI_DUAL_FINAL_SHOW = x with x > 0, Dual Simplex shows the condition + (arcs and nodes) of the network every x iterations. */ + +#define UNIPI_DUAL_FINAL_SHOW 0 + +/* If UNIPI_DUAL_FINAL_SHOW = 1, Dual Simplex shows the final condition + (arcs and nodes) of the network. */ + +#define UNIPI_VIS_DUMMY_ARCS 1 + +/* If UNIPI_VIS_DUMMY_ARCS = 1, Primal Simplex or Dual Simplex shows the + conditions of the dummy arcs. */ + +#define UNIPI_VIS_ARC_UPPER 1 +#define UNIPI_VIS_ARC_COST 1 +#define UNIPI_VIS_ARC_Q_COST 1 +#define UNIPI_VIS_ARC_REDUCT_COST 1 +#define UNIPI_VIS_ARC_STATE 1 +#define UNIPI_VIS_NODE_BASIC_ARC 1 + +/* When Primal Simplex or Dual Simplex shows the conditions of the network, + for every arcs the algorithm shows the flow, for every nodes it shows + balance and potential. These 6 flags decide if the algorithm shows a + particular value of the arcs/nodes; for example if + UNIPI_VIS_ARC_UPPER == 1, the algorithm shows the upper bounds of all + arcs. */ + +#define FOSHOW 0 + +/* If FOSHOW is 1, the algorithm shows the f.o. value every x iterations + (x = UNIPI_PRIMAL_ITER_SHOW or x = UNIPI_DUAL_ITER_SHOW). */ + +#define OPTQUADRATIC 0 + +/* If OPTQUADRATIC is 1 the Primal Simplex, in the quadratic case, tries to + optimize the update of the potential. + Unfortunately this doesn't work well: for this reason it is set to 0. */ + +#ifndef throw +// QUICK HACK +#define throw(x) assert(false); +#endif + +/*--------------------------------------------------------------------------*/ +/*--------------------------- FUNCTIONS ------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +template<class T> +inline T ABS( const T x ) +{ + return( x >= T( 0 ) ? x : - x ); + } + +/*--------------------------------------------------------------------------*/ + +template<class T> +inline void Swap( T &v1 , T &v2 ) +{ + T temp = v1; + v1 = v2; + v2 = temp; + } + +/*--------------------------------------------------------------------------*/ +/*--------------------------- CONSTANTS ------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( QUADRATICCOST == 0 ) + static const char DELETED = -2; // ident for deleted arcs + static const char CLOSED = -1; // ident for closed arcs + static const char BASIC = 0; // ident for basis arcs + static const char AT_LOWER = 1; // ident for arcs in L + static const char AT_UPPER = 2; // ident for arcs in U +#endif + +/* These macros will be used by method MemAllocCandidateList() to set the + values of numCandidateList and hotListSize. There are different macros, + according to: + + - the used Simplex + - the size of the network + - (obviously) the different variables + + This set of values tries to improve the performance of the two algorithms + according to diversified situations. */ + +static const int PRIMAL_LOW_NUM_CANDIDATE_LIST = 30; +static const int PRIMAL_MEDIUM_NUM_CANDIDATE_LIST = 50; +static const int PRIMAL_HIGH_NUM_CANDIDATE_LIST = 200; +static const int PRIMAL_LOW_HOT_LIST_SIZE = 5; +static const int PRIMAL_MEDIUM_HOT_LIST_SIZE = 10; +static const int PRIMAL_HIGH_HOT_LIST_SIZE = 20; +static const int DUAL_LOW_NUM_CANDIDATE_LIST = 6; +static const int DUAL_HIGH_NUM_CANDIDATE_LIST = 10; +static const int DUAL_LOW_HOT_LIST_SIZE = 1; +static const int DUAL_HIGH_HOT_LIST_SIZE = 2; + + + +/*--------------------------------------------------------------------------*/ +/*--------------------------- COSTRUCTOR -----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +MCFSimplex::MCFSimplex( cIndex nmx , cIndex mmx ) + : + MCFClass( nmx , mmx ) +{ + #if( QUADRATICCOST ) + if( numeric_limits<FNumber>::is_integer ) + throw( MCFException( "FNumber must be float if QUADRATICCOST == 1" ) ); + + if( numeric_limits<CNumber>::is_integer ) + throw( MCFException( "CNumber must be float if QUADRATICCOST == 1" ) ); + #endif + + usePrimalSimplex = true; + + newSession = true; + if( nmax && mmax ) + MemAlloc(); + else + nmax = mmax = 0; + + #if( QUADRATICCOST ) + recomputeFOLimits = 100; + // recomputeFOLimits represents the limit of the iteration in which + // quadratic Primal Simplex computes "manually" the f.o. value + EpsOpt = 1e-13; + // EpsOpt is the fixed precision of the quadratic Primal Simplex + #endif + + pricingRule = kCandidateListPivot; + forcedNumCandidateList = 0; + forcedHotListSize = 0; + usePrimalSimplex = true; + nodesP = NULL; + nodesD = NULL; + arcsP = NULL; + arcsD = NULL; + candP = NULL; + candD = NULL; + + if( numeric_limits<CNumber>::is_integer ) + MAX_ART_COST = CNumber( 1e7 ); + else + MAX_ART_COST = CNumber( 1e10 ); + + } // end( MCFSimplex ) + +/*--------------------------------------------------------------------------*/ +/*-------------------------- OTHER INITIALIZATIONS -------------------------*/ +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::LoadNet( cIndex nmx , cIndex mmx , cIndex pn , cIndex pm , + cFRow pU , cCRow pC , cFRow pDfct , + cIndex_Set pSn , cIndex_Set pEn ) +{ + MemDeAllocCandidateList(); + + if( ( nmx != nmax ) || ( mmx != mmax ) ) { + // if the size of the allocated memory changes + + if( nmax && mmax ) { // if the memory was already allocated + MemDeAlloc(true); // deallocate the Primal + MemDeAlloc(false); // and the Dual data structures + nmax = mmax = 0; + } + + if( nmx && mmx ) { // if the new dimensions of the memory are correct + nmax = nmx; + mmax = mmx; + MemAlloc(); + } + } + + if( ( ! nmax ) || ( ! mmax ) ) + // if one of the two dimension of the memory isn't correct + nmax = mmax = 0; + + if( nmax ) { // if the new dimensions of the memory are correct + n = pn; + m = pm; + + if( usePrimalSimplex ) { + stopNodesP = nodesP + n; + dummyRootP = nodesP + nmax; + for( nodePType *node = nodesP ; node != stopNodesP ; node++ ) + node->balance = pDfct[ node - nodesP ]; // initialize nodes + + stopArcsP = arcsP + m; + dummyArcsP = arcsP + mmax; + stopDummyP = dummyArcsP + n; + for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) { + // initialize real arcs + if (pC) + arc->cost = pC[ arc - arcsP ]; + else + arc->cost = 0; + #if( QUADRATICCOST ) + arc->quadraticCost = 0; + #endif + if (pU) + arc->upper = pU[ arc - arcsP ]; + else + arc->upper = Inf<FNumber>(); + arc->tail = nodesP + pSn[ arc - arcsP ] - 1; + arc->head = nodesP + pEn[ arc - arcsP ] - 1; + } + } + else { + stopNodesD = nodesD + n; + dummyRootD = nodesD + nmax; + for( nodeDType *node = nodesD ; node != stopNodesD ; node++ ) + node->balance = pDfct[ node - nodesD ]; // initialize nodes + + stopArcsD = arcsD + m; + dummyArcsD = arcsD + mmax; + stopDummyD = dummyArcsD + n; + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) { + // initialize real arcs + arc->cost = pC[ arc - arcsD ]; + #if( QUADRATICCOST ) + arc->quadraticCost = 0; + #endif + arc->upper = pU[ arc - arcsD ]; + arc->tail = nodesD + pSn[ arc - arcsD ] - 1; + arc->head = nodesD + pEn[ arc - arcsD ] - 1; + } + + CreateAdditionalDualStructures(); + } + + if( pricingRule == kCandidateListPivot ) + MemAllocCandidateList(); + + status = kUnSolved; + } + } // end( MCFSimplex::LoadNet ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::SetAlg( bool UsPrml , char WhchPrc ) +{ + bool newUsePrimalSimplex = UsPrml; + bool oldUsePrimalSimplex = usePrimalSimplex; + char newPricingRule = WhchPrc; + char oldPricingRule = pricingRule; + + usePrimalSimplex = newUsePrimalSimplex; + pricingRule = newPricingRule; + + if( ( ! usePrimalSimplex ) && ( pricingRule == kDantzig) ) { + pricingRule = kFirstEligibleArc; + newPricingRule = pricingRule; + } + + if( ( newUsePrimalSimplex != oldUsePrimalSimplex ) || + ( newPricingRule != oldPricingRule ) ) { + MemDeAllocCandidateList(); + + if( newUsePrimalSimplex != oldUsePrimalSimplex ) { + #if( QUADRATICCOST ) + throw( + MCFException( "Primal Simplex is the only option if QUADRATICCOST == 1" + ) ); + } + #else + MemAlloc(); + nodePType *nP = nodesP; + nodeDType *nD = nodesD; + arcPType *aP = arcsP; + arcDType *aD = arcsD; + + if( newUsePrimalSimplex ) { // from Dual to Primal + if( nodesD == NULL ) + return; + + if( newSession ) + CreateInitialDualBase(); + + dummyRootP = nodesP + nmax; + stopNodesP = nodesP + n; + dummyArcsP = arcsP + mmax; + stopArcsP = arcsP + m; + stopDummyP = dummyArcsP + n; + // Copy the old Dual data structure in a new Primal data structure + while( nD != stopNodesD ) { + nP->prevInT = nodesP + ( nD->prevInT - nodesD ); + nP->nextInT = nodesP + ( nD->nextInT - nodesD ); + nP->enteringTArc = arcsP + ( nD->enteringTArc - arcsD ); + nP->balance = nD->balance; + nP->potential = nD->potential; + nP->subTreeLevel = nD->subTreeLevel; + nP++; + nD++; + } + + dummyRootP->prevInT = NULL; + dummyRootP->nextInT = nodesP + ( dummyRootD->nextInT - nodesD ); + dummyRootP->enteringTArc = arcsP + ( dummyRootD->enteringTArc - arcsD ); + dummyRootP->balance = dummyRootD->balance; + dummyRootP->potential = dummyRootD->potential; + dummyRootP->subTreeLevel = dummyRootD->subTreeLevel; + while( aD != stopArcsD ) { + aP->tail = nodesP + ( aD->tail - nodesD ); + aP->head = nodesP + ( aD->head - nodesD ); + aP->flow = aD->flow; + aP->cost = aD->cost; + aP->ident = aD->ident; + aP->upper = aD->upper; + aP++; + aD++; + } + + aP = dummyArcsP; + aD = dummyArcsD; + while( aD != stopDummyD ) { + aP->tail = nodesP + ( aD->tail - nodesD ); + aP->head = nodesP + ( aD->head - nodesD ); + aP->flow = aD->flow; + aP->cost = aD->cost; + aP->ident = aD->ident; + if( ( ETZ(aP->flow, EpsFlw) ) && (aP->ident == AT_UPPER) ) + aP->ident = AT_LOWER; + + aP->upper = Inf<FNumber>(); + aP++; + aD++; + } + + MemDeAlloc(false); + if( Senstv && ( status != kUnSolved ) ) { + nodePType *node = dummyRootP; + for( int i = 0 ; i < n ; i++ ) + node = node->nextInT; + + node->nextInT = NULL; + dummyRootP->prevInT = NULL; + dummyRootP->enteringTArc = NULL; + // balance the flow + CreateInitialPModifiedBalanceVector(); + PostPVisit( dummyRootP ); + // restore the primal admissibility + BalanceFlow( dummyRootP ); + ComputePotential( dummyRootP ); + } + else + status = kUnSolved; + } + else { // from Primal to Dual + if( nodesP == NULL ) + return; + + if( newSession ) + CreateInitialPrimalBase(); + + dummyRootD = nodesD + nmax; + stopNodesD = nodesD + n; + dummyArcsD = arcsD + mmax; + stopArcsD = arcsD + m; + stopDummyD = dummyArcsD + n; + // Copy the old Primal data structure in a new Dual data structure + while( nP != stopNodesP ) { + nD->prevInT = nodesD + ( nP->prevInT - nodesP ); + nD->nextInT = nodesD + ( nP->nextInT - nodesP ); + nD->enteringTArc = arcsD + ( nP->enteringTArc - arcsP ); + nD->balance = nP->balance; + nD->potential = nP->potential; + nD->subTreeLevel = nP->subTreeLevel; + nP++; + nD++; + } + + dummyRootD->prevInT = NULL; + dummyRootD->nextInT = nodesD + ( dummyRootP->nextInT - nodesP ); + dummyRootD->enteringTArc = NULL; + dummyRootD->balance = dummyRootP->balance; + dummyRootD->potential = dummyRootP->potential; + dummyRootD->subTreeLevel = dummyRootP->subTreeLevel; + while( aP != stopArcsP ) { + aD->tail = nodesD + ( aP->tail - nodesP ); + aD->head = nodesD + ( aP->head - nodesP ); + aD->flow = aP->flow; + aD->cost = aP->cost; + aD->ident = aP->ident; + aD->upper = aP->upper; + aP++; + aD++; + } + + aP = dummyArcsP; + aD = dummyArcsD; + while( aP != stopDummyP ) { + aD->tail = nodesD + ( aP->tail - nodesP ); + aD->head = nodesD + ( aP->head - nodesP ); + aD->flow = aP->flow; + aD->cost = aP->cost; + aD->ident = aP->ident; + aD->upper = 0; + aP++; + aD++; + } + + CreateAdditionalDualStructures(); + MemDeAlloc(true); + nodeDType *node = dummyRootD; + for( int i = 0 ; i < n ; i++ ) + node = node->nextInT; + + node->nextInT = NULL; + dummyRootD->enteringTArc = NULL; + dummyRootD->prevInT = NULL; + if( Senstv && ( status != kUnSolved ) ) { + // fix every flow arc according to its reduct cost + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) { + if( (arc->ident == AT_LOWER) && LTZ( ReductCost( arc ) , EpsCst ) ) { + arc->flow = arc->upper; + arc->ident = AT_UPPER; + } + + if( ( arc->ident == AT_UPPER ) && GTZ( ReductCost( arc ) , EpsCst ) ) { + arc->flow = 0; + arc->ident = AT_LOWER; + } + } + + // balance the flow + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + } + else + status = kUnSolved; + } + //#endif + } +#endif + if( newPricingRule == kFirstEligibleArc ) + if( newUsePrimalSimplex ) + arcToStartP = arcsP; + else + arcToStartD = arcsD; + + if( ( nmax && mmax ) && ( newPricingRule == kCandidateListPivot ) ) + MemAllocCandidateList(); + } + } // end( SetAlg ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::SetPar( int par, int val ) +{ + switch( par ) { + case kAlgPrimal: + if( val == kYes ) + SetAlg( true , pricingRule); + + if( val == kNo ) + SetAlg( false, pricingRule ); + + break; + + case kAlgPricing: + + if( ( val == kDantzig ) || ( val == kFirstEligibleArc ) || + ( val == kCandidateListPivot ) ) + SetAlg( usePrimalSimplex , val ); + + break; + + case kNumCandList: + + MemDeAllocCandidateList(); + forcedNumCandidateList = val; + MemAllocCandidateList(); + forcedNumCandidateList = 0; + forcedHotListSize = 0; + break; + + case kHotListSize: + + MemDeAllocCandidateList(); + forcedHotListSize = val; + MemAllocCandidateList(); + forcedNumCandidateList = 0; + forcedHotListSize = 0; + break; + + case kRecomputeFOLimits: + + recomputeFOLimits = val; + break; + + default: + + MCFClass::SetPar(par, val); + } + } // end( SetPar( int ) ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::SetPar( int par , double val ) +{ + switch( par ) { + case kEpsOpt: + + EpsOpt = val; + break; + + default: + MCFClass::SetPar( par , val ); + } + } // end( SetPar( double ) + +/*-------------------------------------------------------------------------*/ +/*--------------- METHODS FOR SOLVING THE PROBLEM -------------------------*/ +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::SolveMCF( void ) +{ + if( MCFt ) + MCFt->Start(); + + //if( status == kUnSolved ) + if(usePrimalSimplex ) + CreateInitialPrimalBase(); + else + CreateInitialDualBase(); + + newSession = false; + if( usePrimalSimplex ) + PrimalSimplex(); + else + DualSimplex(); + + if( MCFt ) + MCFt->Stop(); + + } // end( MCFSimplex::SolveMCF ) + +/*--------------------------------------------------------------------------*/ +/*---------------------- METHODS FOR READING RESULTS -----------------------*/ +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MCFGetX( FRow F , Index_Set nms , cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + if( usePrimalSimplex ) + for( Index i = strt ; i < stp ; i++ ) { + FNumber tXi = ( arcsP + i )->flow; + if( GTZ( tXi , EpsFlw ) ) { + *(F++) = tXi; + *(nms++) = i; + } + } + else + for( Index i = strt ; i < stp ; i++ ) { + FNumber tXi = ( arcsD + i )->flow; + if( GTZ( tXi , EpsFlw ) ) { + *(F++) = tXi; + *(nms++) = i; + } + } + + *nms = Inf<Index>(); + } + else + if( usePrimalSimplex ) + for( Index i = strt; i < stp; i++ ) + *(F++) = ( arcsP + i )->flow; + else + for( Index i = strt; i < stp; i++ ) + *(F++) = ( arcsD + i )->flow; + + } // end( MCFSimplex::MCFGetX( some ) ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MCFGetRC( CRow CR , cIndex_Set nms , cIndex strt , Index stp ) +{ + if( nms ) { + while( *nms < strt ) + nms++; + + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(CR++) = CNumber( ReductCost( arcsP + h ) ); + else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(CR++) = ReductCost( arcsD + h ); + } + else { + if( stp > m ) + stp = m; + + if( usePrimalSimplex ) + for( arcPType* arc = arcsP + strt ; arc < arcsP + stp ; arc++ ) + *(CR++) = CNumber( ReductCost( arc ) ); + else + for( arcDType* arc = arcsD + strt ; arc < arcsD + stp ; arc++ ) + *(CR++) = ReductCost( arc ); + } + } // end( MCFSimplex::MCFGetRC( some ) ) + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::CNumber MCFSimplex::MCFGetRC( cIndex i ) +{ + if( usePrimalSimplex ) + return CNumber( ReductCost( arcsP + i ) ); + else + return( ReductCost( arcsD + i ) ); + + } // end( MCFSimplex::MCFGetRC( i ) ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MCFGetPi( CRow P , cIndex_Set nms , cIndex strt , Index stp ) +{ + if( stp > n ) + stp = n; + + if( nms ) { + while( *nms < strt ) + nms++; + + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(P++) = CNumber( (nodesP + h)->potential ); + else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(P++) = (nodesD + h )->potential; + } + else + if( usePrimalSimplex ) + for( nodePType *node = nodesP + strt ; node < ( nodesP + stp ) ; node++ ) + *(P++) = CNumber( node->potential ); + else + for( nodeDType *node = nodesD + strt ; node++ < ( nodesD + stp ) ; node++ ) + *(P++) = node->potential; + + } // end( MCFSimplex::MCFGetPi( some ) ) + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::FONumber MCFSimplex::MCFGetFO( void ) +{ + if( status == kOK ) + return( (FONumber) GetFO() ); + else + if( status == kUnbounded ) + return( - Inf<FONumber>() ); + else + return( Inf<FONumber>() ); + + } // end( MCFSimplex::MCFGetFO ) + +/*-------------------------------------------------------------------------*/ +/*----------METHODS FOR READING THE DATA OF THE PROBLEM--------------------*/ +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::MCFArcs( Index_Set Startv , Index_Set Endv , + cIndex_Set nms , cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) + nms++; + + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) { + if( Startv ) + *(Startv++) = Index( (arcsP + h)->tail - nodesP) - USENAME0; + if( Endv ) + *(Endv++) = Index( (arcsP + h)->head - nodesP ) - USENAME0; + } + else + for( Index h ; ( h = *(nms++) ) < stp ; ) { + if( Startv ) + *(Startv++) = Index( (arcsD + h)->tail - nodesD) - USENAME0; + if( Endv ) + *(Endv++) = Index( (arcsD + h)->head - nodesD ) - USENAME0; + } + } + else + if( usePrimalSimplex ) + for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ ) { + if( Startv ) + *(Startv++) = Index( arc->tail - nodesP ) - USENAME0; + if( Endv ) + *(Endv++) = Index( arc->head - nodesP ) - USENAME0; + } + else + for( arcDType* arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ ) { + if( Startv ) + *(Startv++) = Index( arc->tail - nodesD ) - USENAME0; + if( Endv ) + *(Endv++) = Index( arc->head - nodesD ) - USENAME0; + } + + } // end( MCFSimplex::MCFArcs ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::MCFCosts( CRow Costv , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) + nms++; + + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Costv++) = (arcsP + h)->cost; + else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Costv++) = (arcsD + h)->cost; + } + else + if( usePrimalSimplex ) + for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ ) + *(Costv++) = arc->cost; + else + for( arcDType* arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ ) + *(Costv++) = arc->cost; + + } // end( MCFSimplex::MCFCosts ( some ) ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::MCFQCoef( CRow Qv , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) + nms++; + + #if( QUADRATICCOST ) + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Qv++) = (arcsP + h)->quadraticCost; + else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Qv++) = (arcsD + h)->quadraticCost; + #else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Qv++) = 0; + #endif + } + else + #if( QUADRATICCOST ) + if( usePrimalSimplex ) + for( arcPType* arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ ) + *(Qv++) = arc->quadraticCost; + else + for( arcDType* arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ ) + *(Qv++) = arc->quadraticCost; + #else + for( Index h = strt ; h++ < stp ; ) + *(Qv++) = 0; + #endif + + } // end( MCFSimplex::MCFQCoef ( some ) ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::MCFUCaps( FRow UCapv , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) + nms++; + + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(UCapv++) = (arcsP + h)->upper; + else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(UCapv++) = (arcsD + h)->upper; + } + else + if( usePrimalSimplex ) + for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp ) ; arc++ ) + *(UCapv++) = arc->upper; + else + for( arcDType* arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ ) + *(UCapv++) = arc->upper; + + } // end( MCFSimplex::MCFUCaps ( some ) ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::MCFDfcts( FRow Dfctv , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > n ) + stp = n; + + if( nms ) { + while( *nms < strt ) + nms++; + + if( usePrimalSimplex ) + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Dfctv++) = ( nodesP + h )->balance; + else + for( Index h ; ( h = *(nms++) ) < stp ; ) + *(Dfctv++) = (nodesD + h )->balance; + } + else + if( usePrimalSimplex ) + for( nodePType* node = nodesP + strt ; node < ( nodesP + stp ) ; node++ ) + *(Dfctv++) = node->balance; + else + for( nodeDType* node = nodesD + strt ; node < ( nodesD + stp ) ; node++ ) + *(Dfctv++) = node->balance; + + } // end( MCFSimplex::MCFDfcts ) + +/*-------------------------------------------------------------------------*/ +/*--------- METHODS FOR ADDING / REMOVING / CHANGING DATA -----------------*/ +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgCosts( cCRow NCost , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) { + nms++; + NCost++; + } + + cIndex_Set tnms = nms; // nms may be needed below + if( usePrimalSimplex ) + for( Index h ; ( h = *(tnms++) ) < stp ; ) + arcsP[ h ].cost = *(NCost++); + else + for( Index h ; ( h = *(tnms++) ) < stp ; ) + arcsD[ h ].cost = *(NCost++); + } + else + if( usePrimalSimplex ) + for( arcPType *arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ ) + arc->cost = *(NCost++); + else + for( arcDType *arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ ) + arc->cost = *(NCost++); + + if( Senstv && ( status != kUnSolved ) ) + if( usePrimalSimplex ) + ComputePotential( dummyRootP ); + else { + #if( QUADRATICCOST == 0 ) + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) + if( arc->ident > BASIC ) + if( GTZ( ReductCost( arc ) , EpsCst ) ) { + arc->flow = 0; + arc->ident = AT_LOWER; + } + else { + arc->flow = arc->upper; + arc->ident = AT_UPPER; + } + + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + #endif + } + else + status = kUnSolved; + + } // end( MCFSimplex::ChgCosts ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgCost( Index arc , cCNumber NCost ) +{ + if( arc >= m ) + return; + + if( usePrimalSimplex ) + ( arcsP + arc )->cost = NCost; + else + ( arcsD + arc )->cost = NCost; + + if( Senstv && ( status != kUnSolved ) ) { + if( usePrimalSimplex ) { + nodePType *node = ( arcsP + arc )->tail; + if( ( ( arcsP + arc )->head)->subTreeLevel < node->subTreeLevel ) + node = ( arcsP + arc )->head; + + ComputePotential( dummyRootP ); + } + else { + #if( QUADRATICCOST == 0 ) + nodeDType *node = ( arcsD + arc )->tail; + if( ( ( arcsD + arc )->head )->subTreeLevel < node->subTreeLevel ) + node = ( arcsD + arc )->head; + + ComputePotential( dummyRootD ); + for( arcDType *a = arcsD ; a != stopArcsD ; a++) + if( a->ident > BASIC ) + if( GTZ( ReductCost( a ) , EpsCst ) ) { + a->flow = 0; + a->ident = AT_LOWER; + } + else { + a->flow = a->upper; + a->ident = AT_UPPER; + } + + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + #endif + } + } + else + status = kUnSolved; + + } // end( MCFSimplex::ChgCost ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgQCoef( cCRow NQCoef , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + #if( QUADRATICCOST ) + if( nms ) { + while( *nms < strt ) { + nms++; + NQCoef++; + } + + cIndex_Set tnms = nms; // nms may be needed below + if( usePrimalSimplex ) + for( Index h ; ( h = *(tnms++) ) < stp ; ) + arcsP[ h ].quadraticCost = *(NQCoef++); + else + for( Index h ; ( h = *(tnms++) ) < stp ; ) + arcsD[ h ].quadraticCost = *(NQCoef++); + } + else + if( usePrimalSimplex ) + for( arcPType *arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ ) + arc->quadraticCost = *(NQCoef++); + else + for( arcDType *arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ ) + arc->quadraticCost = *(NQCoef++); + + if( Senstv && (status != kUnSolved ) ) + ComputePotential( dummyRootP ); + else + status = kUnSolved; + #else + if( NQCoef ) + throw( MCFException( "ChgQCoef: nonzero coefficients not allowed" ) ); + #endif +} // end( MCFSimplex::ChgQCoef ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgQCoef( Index arc , cCNumber NQCoef ) +{ + #if( QUADRATICCOST ) + if( arc >= m ) + return; + + if( usePrimalSimplex ) + ( arcsP + arc )->quadraticCost = NQCoef; + else + ( arcsD + arc )->quadraticCost = NQCoef; + + if( Senstv && ( status != kUnSolved ) ) { + nodePType *node = ( arcsP + arc )->tail; + if( ( ( arcsP + arc )->head )->subTreeLevel < node->subTreeLevel ) + node = ( arcsP + arc )->head; + + ComputePotential( node ); + } + else + status = kUnSolved; + #else + if( NQCoef ) + throw( MCFException( "ChgQCoef: nonzero coefficients not allowed" ) ); + #endif + } // end( MCFSimplex::ChgQCoef ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgDfcts( cFRow NDfct , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) { + nms++; + NDfct++; + } + + cIndex_Set tnms = nms; // nms may be needed below + if( usePrimalSimplex ) + for( Index h ; ( h = *(tnms++) ) < stp ; ) + nodesP[ h ].balance = *(NDfct++); + else + for( Index h ; ( h = *(tnms++) ) < stp ; ) + nodesD[ h ].balance = *(NDfct++); + } + else + if( usePrimalSimplex ) + for( nodePType *node = nodesP + strt ; node < ( nodesP + stp ) ; node++ ) + node->balance = *(NDfct++); + else + for( nodeDType *node = nodesD + strt ; node < ( nodesD + stp ) ; node++ ) + node->balance = *(NDfct++); + + if( Senstv && (status != kUnSolved ) ) + if( usePrimalSimplex ) { + CreateInitialPModifiedBalanceVector(); + PostPVisit( dummyRootP ); + BalanceFlow( dummyRootP ); + ComputePotential( dummyRootP ); + } + else { + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + } + else + status = kUnSolved; + + } // end( MCFSimplex::ChgDfcts ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgDfct( Index nod , cFNumber NDfct ) +{ + if( nod > n ) + return; + + if( usePrimalSimplex ) + ( nodesP + nod - 1 )->balance = NDfct; + else + ( nodesD + nod - 1 )->balance = NDfct; + + if( Senstv && (status != kUnSolved ) ) + if( usePrimalSimplex ) { + CreateInitialPModifiedBalanceVector(); + PostPVisit( dummyRootP ); + BalanceFlow( dummyRootP ); + ComputePotential( dummyRootP ); + } + else { + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + } + else + status = kUnSolved; + + } // end( MCFSimplex::ChgDfct ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgUCaps( cFRow NCap , cIndex_Set nms , + cIndex strt , Index stp ) +{ + if( stp > m ) + stp = m; + + if( nms ) { + while( *nms < strt ) { + nms++; + NCap++; + } + + cIndex_Set tnms = nms; // nms may be needed below + if( usePrimalSimplex ) + for( Index h ; ( h = *(tnms++) ) < stp ; ) + arcsP[ h ].upper = *(NCap++); + else + for( Index h ; ( h = *(tnms++) ) < stp ; ) + arcsD[ h ].upper = *(NCap++); + } + else + if( usePrimalSimplex ) + for( arcPType *arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ ) + arc->upper = *(NCap++); + else + for( arcDType *arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ ) + arc->upper = *(NCap++); + + if( Senstv && (status != kUnSolved ) ) { + if( usePrimalSimplex ) { + for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++) + #if( QUADRATICCOST ) + if( GT( arc->flow , arc->upper , EpsFlw ) ) + arc->flow = arc->upper; + #else + if( GT(arc->flow , arc->upper , EpsFlw ) || + ( ( arc->ident == AT_UPPER ) && + ( ! ETZ( arc->flow - arc->upper , EpsFlw ) ) ) ) + arc->flow = arc->upper; + #endif + + CreateInitialPModifiedBalanceVector(); + PostPVisit( dummyRootP ); + BalanceFlow( dummyRootP ); + ComputePotential( dummyRootP ); + } + else { + #if( QUADRATICCOST == 0 ) + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) + if( ( GT( arc->flow , arc->upper , EpsFlw ) && ( arc->ident != BASIC ) ) || + ( ( arc->ident == AT_UPPER ) && + ( ! ETZ( arc->flow - arc->upper , EpsFlw ) ) ) ) + arc->flow = arc->upper; + + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + ComputePotential( dummyRootD ); + #endif + } + } + else + status = kUnSolved; + + } // end( MCFSimplex::ChgUCaps ) + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::ChgUCap( Index arc , cFNumber NCap ) +{ + if( arc >= m ) + return; + + if( usePrimalSimplex ) + ( arcsP + arc )->upper = NCap; + else + ( arcsD + arc )->upper = NCap; + + if( Senstv && (status != kUnSolved ) ) { + if( usePrimalSimplex ) { + #if( QUADRATICCOST ) + if( GT( ( arcsP + arc )->flow , ( arcsP + arc )->upper , EpsFlw ) ) + ( arcsP + arc )->flow = ( arcsP + arc )->upper; + #else + if( GT( ( arcsP + arc )->flow , ( arcsP + arc )->upper , EpsFlw ) || + ( ( ( arcsP + arc )->ident == AT_UPPER ) && + ( ! ETZ( ( arcsP + arc )->flow - ( arcsP + arc )->upper , EpsFlw ) ) ) ) + ( arcsP + arc )->flow = ( arcsP + arc )->upper; + #endif + + CreateInitialPModifiedBalanceVector(); + PostPVisit( dummyRootP ); + BalanceFlow( dummyRootP ); + ComputePotential( dummyRootP ); + } + else { + #if( QUADRATICCOST == 0 ) + if( ( GT( ( arcsD + arc )->flow , ( arcsD + arc )->upper , EpsFlw ) && + ( ( ( arcsD + arc )->ident != BASIC ) ) ) || + ( ( ( arcsD + arc )->ident == AT_UPPER ) && + ( ! ETZ( ( arcsD + arc )->flow - ( arcsD + arc )->upper , EpsFlw ) ) ) ) { + ( arcsD + arc )->flow = ( arcsD + arc )->upper; + ( arcsD + arc )->ident = AT_UPPER; + } + + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + ComputePotential( dummyRootD ); + #endif + } + } + else + status = kUnSolved; + + } // end( MCFSimplex::ChgUCap ) + +/*-------------------------------------------------------------------------*/ + +bool MCFSimplex::IsClosedArc( cIndex name ) +{ + if( name >= m ) + return( false ); + + #if( QUADRATICCOST ) + return( ( arcsP + name )->cost == Inf<CNumber>() ); + #else + if( usePrimalSimplex ) + return( ( ( arcsP + name )->ident < BASIC) ); + else + return( ( ( arcsD + name )->ident < BASIC) ); + #endif + } + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::CloseArc( cIndex name ) +{ + if( name >= m ) + return; + + if( usePrimalSimplex ) { + arcPType *arc = arcsP + name; + #if( QUADRATICCOST ) + if( arc->cost == Inf<CNumber>() ) + return; + + arc->cost = Inf<CNumber>(); + #else + if( arc->ident < BASIC ) + return; + + arc->ident = CLOSED; + #endif + + arc->flow = 0; + + if( Senstv && ( status != kUnSolved ) ) { + nodePType *node = NULL; + if( (arc->tail)->enteringTArc == arc ) + node = arc->tail; + + if( (arc->head)->enteringTArc == arc ) + node = arc->head; + + if( node ) { + node->enteringTArc = dummyArcsP + ( node - nodesP ); + nodePType *last = CutAndUpdateSubtree( node , -node->subTreeLevel + 1 ); + PasteSubtree( node , last , dummyRootP ); + node->enteringTArc = dummyArcsP + ( node - nodesP ); + } + + CreateInitialPModifiedBalanceVector(); + PostPVisit(dummyRootP); + BalanceFlow(dummyRootP); + ComputePotential(dummyRootP); + } + else + status = kUnSolved; + } + else { + #if( QUADRATICCOST == 0 ) + arcDType *arc = arcsD + name; + if( arc->ident < BASIC ) + return; + + arc->flow = 0; + arc->ident = CLOSED; + + if( Senstv && ( status != kUnSolved ) ) { + nodeDType *node = NULL; + if( ( arc->tail )->enteringTArc == arc) + node = arc->tail; + + if( ( arc->head )->enteringTArc == arc ) + node = arc->head; + + if( node ) { + node->enteringTArc = dummyArcsD + ( node - nodesD ); + nodeDType *last = CutAndUpdateSubtree( node , -node->subTreeLevel + 1 ); + PasteSubtree( node , last , dummyRootD ); + node->enteringTArc = dummyArcsD + ( node - nodesD ); + ComputePotential( dummyRootD ); + + for( arcDType *a = arcsD ; a != stopArcsD ; a++ ) + if( a->ident > BASIC ) + if( GTZ( ReductCost( a ) , EpsCst ) ) { + a->flow = 0; + a->ident = AT_LOWER; + } + else { + a->flow = a->upper; + a->ident = AT_UPPER; + } + } + + CreateInitialDModifiedBalanceVector(); + PostDVisit(dummyRootD); + ComputePotential(dummyRootD); + } + else + status = kUnSolved; + #endif + } + } // end( MCFSimplex::CloseArc ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::DelNode( cIndex name ) +{ + if( name >= n ) + return; + + if( usePrimalSimplex ) { + nodePType *node = nodesP + name; + nodePType *last = CutAndUpdateSubtree(node, -node->subTreeLevel); + nodePType *n = node->nextInT; + while( n ) { + if( n->subTreeLevel == 1 ) + n->enteringTArc = dummyArcsP + ( n - nodesP ); + + n = n->nextInT; + } + + PasteSubtree( node , last , dummyRootP ); + n = node->nextInT; + dummyRootP->nextInT = n; + n->prevInT = dummyRootP; + + for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) { + if( ( ( arc->tail ) == node) || ( ( arc->head ) == node ) ) { + arc->flow = 0; + #if( QUADRATICCOST ) + arc->cost = Inf<CNumber>(); + #else + arc->ident = CLOSED; + #endif + } + } + + for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) { + if( ( ( arc->tail ) == node ) || ( ( arc->head ) == node ) ) { + arc->flow = 0; + #if( QUADRATICCOST ) + arc->cost = Inf<CNumber>(); + #else + arc->ident = CLOSED; + #endif + } + } + + CreateInitialPModifiedBalanceVector(); + PostPVisit( dummyRootP ); + BalanceFlow( dummyRootP ); + ComputePotential( dummyRootP ); + } + else { + #if( QUADRATICCOST == 0 ) + nodeDType *node = nodesD + name; + nodeDType *last = CutAndUpdateSubtree( node , -node->subTreeLevel ); + nodeDType *n = node->nextInT; + while( n ) { + if( n->subTreeLevel == 1 ) + n->enteringTArc = dummyArcsD + ( n - nodesD ); + + n = n->nextInT; + } + + PasteSubtree( node , last , dummyRootD ); + n = node->nextInT; + dummyRootD->nextInT = n; + n->prevInT = dummyRootD; + + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) + if( ( ( arc->tail ) == node) || ( ( arc->head ) == node ) ) { + arc->flow = 0; + arc->ident = CLOSED; + } + + for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ ) + if( ( ( arc->tail ) == node ) || ( ( arc->head ) == node ) ) { + arc->flow = 0; + arc->ident = CLOSED; + } + + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + ComputePotential( dummyRootD ); + #endif + } + } // end( MCFSimplex::DelNode ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::OpenArc( cIndex name ) +{ + if( name >= m ) + return; + + if( usePrimalSimplex ) { + /* Quadratic case is not implemented for a theory bug. + Infact a closed arc in the quadratic case has its cost fixed to infinity, and + it's impossible to restore the old value. */ + + #if( QUADRATICCOST == 0 ) + arcPType *arc = arcsP + name; + if( arc->ident == CLOSED ) { + arc->ident = AT_LOWER; + arc->flow = 0; + } + #endif + } + else { + #if( QUADRATICCOST == 0 ) + arcDType *arc = arcsD + name; + if( arc->ident == CLOSED ) { + if( GTZ( ReductCost( arc ) , EpsCst ) ) + arc->ident = AT_LOWER; + else { + arc->ident = AT_UPPER; + arc->flow = arc->upper; + if( Senstv && ( status != kUnSolved ) ) { + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + } + else + status = kUnSolved; + } + } + #endif + } + } // end( MCFSimplex:OpenArc ) + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::Index MCFSimplex::AddNode( cFNumber aDfct ) +{ + if( n >= nmax ) + return( Inf<Index>() ); + + n++; + if( usePrimalSimplex ) { + nodePType *newNode = nodesP + n - 1; + stopArcsP->tail = newNode; + stopArcsP->head = dummyRootP; + stopArcsP->upper = Inf<FNumber>(); + stopArcsP->flow = 0; + stopArcsP->cost = Inf<CNumber>(); + #if( QUADRATICCOST ) + stopArcsP->quadraticCost = 0; + #else + stopArcsP->ident = BASIC; + #endif + stopArcsP++; + newNode->balance = aDfct; + newNode->prevInT = dummyRootP; + newNode->nextInT = dummyRootP->nextInT; + (dummyRootP->nextInT)->prevInT = newNode; + dummyRootP->nextInT = newNode; + newNode->enteringTArc = stopArcsP--; + newNode->potential = 0; + #if(QUADRATICCOST) + newNode->sumQuadratic = 0; + #endif + } + else { + #if( QUADRATICCOST == 0 ) + nodeDType *newNode = nodesD + n - 1; + stopArcsD->tail = newNode; + stopArcsD->head = dummyRootD; + stopArcsD->upper = 0; + stopArcsD->flow = 0; + stopArcsD->cost = Inf<CNumber>(); + stopArcsD->ident = BASIC; + newNode->balance = aDfct; + newNode->prevInT = dummyRootD; + newNode->nextInT = dummyRootD->nextInT; + (dummyRootD->nextInT)->prevInT = newNode; + dummyRootD->nextInT = newNode; + newNode->enteringTArc = stopArcsD; + newNode->potential = 0; + newNode->firstFs = stopArcsD; + newNode->firstBs = NULL; + stopArcsD->nextFs = NULL; + stopArcsD->nextBs = dummyRootD->firstBs; + dummyRootD->firstBs = stopArcsD; + stopArcsD++; + #endif + } + + return( n ); + + } // end( MCFSimplex::AddNode ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::ChangeArc( cIndex name , cIndex nSN , cIndex nEN ) +{ + if( name >= m ) + return; + + CloseArc( name ); + + if( usePrimalSimplex ) { + if( nSN <= n ) + (arcsP + name)->tail = (nodesP + nSN + USENAME0 - 1); + if( nEN <= n ) + (arcsP + name)->head = (nodesP + nEN + USENAME0 - 1); + } + else { + if( nSN <= n ) + (arcsD + name)->tail = (nodesD + nSN + USENAME0 - 1); + if( nEN <= n ) + (arcsD + name)->head = (nodesD + nEN + USENAME0 - 1); + } + + OpenArc( name ); + + } // end( MCFSimplex::ChangeArc ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::DelArc( cIndex name ) +{ + if( name >= m ) + return; + + if( usePrimalSimplex ) { + arcPType *arc = arcsP + name; + #if( QUADRATICCOST ) + if( arc->upper == -Inf<FNumber>() ) + return; + + if( arc->cost < Inf<CNumber>() ) + #else + if( arc->cost == DELETED ) + return; + + if( arc->cost >= BASIC ) + #endif + CloseArc( name ); + + #if( QUADRATICCOST ) + arc->upper = -Inf<FNumber>(); + #else + arc->ident = DELETED; + #endif + } + else { + #if( QUADRATICCOST == 0 ) + arcDType *arc = arcsD + name; + if( arc->cost == DELETED ) + return; + + if( arc->cost >= BASIC ) + CloseArc( name ); + + arc->ident = DELETED; + #endif + } + } // end( MCFSimplex::DelArc ) + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::Index MCFSimplex::AddArc( cIndex Start , cIndex End , + cFNumber aU , cCNumber aC ) +{ + if( usePrimalSimplex ) { + arcPType *arc = arcsP; + #if( QUADRATICCOST ) + while( ( arc->upper != -Inf<FNumber>() ) && ( arc != stopArcsP ) ) + #else + while( ( arc->ident > DELETED ) && ( arc != stopArcsP ) ) + #endif + arc++; + + if( arc == stopArcsP ) { + if( m >= mmax ) + return( Inf<Index>() ); + + m++; + stopArcsP++; + } + + Index pos = ( arc - arcsP ) + 1; + arc->tail = nodesP + Start + USENAME0 - 1; + arc->head = nodesP + End + USENAME0 - 1; + arc->upper = aU; + arc->cost = aC; + arc->flow = 0; + #if( QUADRATICCOST ) + arc->quadraticCost = 0; + #else + arc->ident = AT_LOWER; + #endif + ComputePotential( dummyRootP ); + return( pos ); + } + else { + #if( QUADRATICCOST == 0 ) + arcDType *arc = arcsD; + while( ( arc->ident > DELETED ) && ( arc != stopArcsD ) ) + arc++; + + if( arc == stopArcsD ) { + if( m >= mmax ) + return( Inf<Index>() ); + + m++; + stopArcsD++; + } + + Index pos = ( arc - arcsD ) + 1; + arc->tail = nodesD + Start + USENAME0 - 1; + arc->head = nodesD + End + USENAME0 - 1; + arc->upper = aU; + arc->cost = aC; + if( GEZ( ReductCost( arc ) , EpsCst ) ) { + arc->flow = 0; + arc->ident = AT_LOWER; + if( Senstv && ( status != kUnSolved ) ) { + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + ComputePotential( dummyRootD ); + } + else + status = kUnSolved; + } + else { + arc->flow = arc->upper; + arc->ident = AT_UPPER; + if( Senstv && ( status != kUnSolved ) ) { + CreateInitialDModifiedBalanceVector(); + PostDVisit( dummyRootD ); + ComputePotential( dummyRootD ); + } + else + status = kUnSolved; + } + + ComputePotential( dummyRootD ); + return( pos ); + #endif + } + } // end( MCFSimplex::AddArc ) + +/*--------------------------------------------------------------------------*/ + +bool MCFSimplex::IsDeletedArc( cIndex name ) +{ + if( name >= m ) + return( false ); + + #if( QUADRATICCOST ) + return( ( ( arcsP + name )->upper == -Inf<FNumber>() ) ); + #else + if( usePrimalSimplex ) + return( ( arcsP + name )->ident == DELETED ); + else + return( ( arcsD + name )->ident == DELETED ); + #endif + } + +/*--------------------------------------------------------------------------*/ +/*------------------------------ DESTRUCTOR --------------------------------*/ +/*--------------------------------------------------------------------------*/ + +MCFSimplex::~MCFSimplex() +{ + MemDeAllocCandidateList(); + MemDeAlloc(true); + MemDeAlloc(false); + } + +/*--------------------------------------------------------------------------*/ +/*---------------------------- PRIVATE METHODS -----------------------------*/ +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MemAlloc( void ) +{ + if( usePrimalSimplex ) { + nodesP = new nodePType[ nmax + 1 ]; // array of nodes + arcsP = new arcPType[ mmax + nmax ]; // array of arcs + dummyArcsP = arcsP + mmax; // artificial arcs are in the last + // nmax positions of the array arcs[] + } + else { + nodesD = new nodeDType[ nmax + 1 ]; // array of nodes + arcsD = new arcDType[ mmax + nmax ]; // array of arcs + dummyArcsD = arcsD + mmax; // artificial arcs are in the last nmax + // positions of the array arcs[] + } + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MemDeAlloc( bool whatDeAlloc ) +{ + if( whatDeAlloc ) { + delete[] nodesP; + delete[] arcsP; + nodesP = NULL; + arcsP = NULL; + } + else { + delete[] nodesD; + delete[] arcsD; + nodesD = NULL; + arcsD = NULL; + } + MemDeAllocCandidateList( ); + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MemAllocCandidateList( void ) +{ + if( usePrimalSimplex ) { + if( m < 10000 ) { + numCandidateList = PRIMAL_LOW_NUM_CANDIDATE_LIST; + hotListSize = PRIMAL_LOW_HOT_LIST_SIZE; + } + else + if( m > 100000 ) { + numCandidateList = PRIMAL_HIGH_NUM_CANDIDATE_LIST; + hotListSize = PRIMAL_HIGH_HOT_LIST_SIZE ; + } + else { + numCandidateList = PRIMAL_MEDIUM_NUM_CANDIDATE_LIST; + hotListSize = PRIMAL_MEDIUM_HOT_LIST_SIZE; + } + + #if( QUADRATICCOST ) + int coef = 1; + // If the number of the arcs is more than 10000, numCandidateList and hotListSize + // are increased to improve the performance of the Quadratic Primal Simplex + if( m > 10000 ) + coef = 10; + + numCandidateList = numCandidateList * coef; + hotListSize = hotListSize * coef; + #endif + + if( forcedNumCandidateList > 0 ) + numCandidateList = forcedNumCandidateList; + + if( forcedHotListSize > 0 ) + hotListSize = forcedHotListSize; + + candP = new primalCandidType[ hotListSize + numCandidateList + 1 ]; + } + else { + if( n < 10000 ) { + numCandidateList = DUAL_LOW_NUM_CANDIDATE_LIST; + hotListSize = DUAL_LOW_HOT_LIST_SIZE; + } + else { + numCandidateList = DUAL_HIGH_NUM_CANDIDATE_LIST; + hotListSize = DUAL_HIGH_HOT_LIST_SIZE; + } + + if( forcedNumCandidateList > 0 ) + numCandidateList = forcedNumCandidateList; + + if( forcedHotListSize > 0 ) + hotListSize = forcedHotListSize; + + candD = new dualCandidType[ hotListSize + numCandidateList + 1 ]; + } + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::MemDeAllocCandidateList( void ) +{ + delete[] candP; + candP = NULL; + delete[] candD; + candD = NULL; +} + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::CreateInitialPrimalBase( void ) +{ + arcPType *arc; + nodePType *node; + for( arc = arcsP ; arc != stopArcsP ; arc++ ) { // initialize real arcs + arc->flow = 0; + #if( QUADRATICCOST == 0 ) + arc->ident = AT_LOWER; + #endif + } + + for( arc = dummyArcsP ; arc != stopDummyP ; arc++ ) { // initialize dummy arcs + node = nodesP + ( arc - dummyArcsP ); + if( node->balance > 0 ) { // sink nodes + arc->tail = dummyRootP; + arc->head = node; + arc->flow = node->balance; + } + else { // source nodes or transit node + arc->tail = node; + arc->head = dummyRootP; + arc->flow = -node->balance; + } + + arc->cost = MAX_ART_COST; + #if( QUADRATICCOST ) + arc->quadraticCost = 0; + #else + arc->ident = BASIC; + #endif + arc->upper = Inf<FNumber>(); + } + + dummyRootP->balance = 0; + dummyRootP->prevInT = NULL; + dummyRootP->nextInT = nodesP; + dummyRootP->enteringTArc = NULL; + #if( QUADRATICCOST ) + dummyRootP->sumQuadratic = 0; + #endif + dummyRootP->potential = MAX_ART_COST; + dummyRootP->subTreeLevel = 0; + for( node = nodesP ; node != stopNodesP ; node++) { // initialize nodes + node->prevInT = node - 1; + node->nextInT = node + 1; + node->enteringTArc = dummyArcsP + (node - nodesP); + #if( QUADRATICCOST ) + node->sumQuadratic = (node->enteringTArc)->quadraticCost; + #endif + if( node->balance > 0 ) // sink nodes + node->potential = 2 * MAX_ART_COST; + else // source nodes or transit node + node->potential = 0; + + node->subTreeLevel = 1; + } + + nodesP->prevInT = dummyRootP; + ( nodesP + n - 1 )->nextInT = NULL; + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::CreateInitialDualBase( void ) +{ + arcDType *arc; + nodeDType *node; + for( arc = dummyArcsD ; arc != stopDummyD ; arc++ ) { // initialize dummy arcs + node = nodesD + ( arc - dummyArcsD ); + arc->tail = node; + arc->head = dummyRootD; + arc->flow = -node->balance; + arc->cost = MAX_ART_COST; + #if( QUADRATICCOST ) + arc->quadraticCost = 0; + #else + arc->ident = BASIC; + #endif + arc->upper = 0; + } + + for( arc = arcsD ; arc != stopArcsD ; arc++ ) { // initialize real arcs + if( GTZ( arc->cost , EpsCst ) ) { + arc->flow = 0; + #if( QUADRATICCOST == 0 ) + arc->ident = AT_LOWER; + #endif + } + else { + #if( QUADRATICCOST == 0 ) + arc->ident = AT_UPPER; + #endif + arc->flow = arc->upper; + ( dummyArcsD + ( ( arc->tail ) - nodesD ) )->flow = + ( dummyArcsD + ( ( arc->tail ) - nodesD ) )->flow - arc->upper; + + ( dummyArcsD + ( ( arc->head ) - nodesD ) )->flow = + ( dummyArcsD + ( ( arc->head ) - nodesD ) )->flow + arc->upper; + } + } + + dummyRootD->balance = 0; + dummyRootD->prevInT = NULL; + dummyRootD->nextInT = nodesD; + dummyRootD->enteringTArc = NULL; + #if( QUADRATICCOST ) + dummyRootD->sumQuadratic = 0; + #endif + dummyRootD->potential = MAX_ART_COST; + dummyRootD->subTreeLevel = 0; + for( node = nodesD ; node != stopNodesD ; node++ ) { // initialize nodes + node->prevInT = node - 1; + node->nextInT = node + 1; + node->enteringTArc = dummyArcsD + ( node - nodesD ); + #if( QUADRATICCOST ) + node->sumQuadratic = ( node->enteringTArc )->quadraticCost; + #endif + node->potential = 0; + node->subTreeLevel = 1; + node->whenInT2 = 0; + } + + nodesD->prevInT = dummyRootD; + ( nodesD + n - 1 )->nextInT = NULL; + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::CreateAdditionalDualStructures( void ) +{ + // this method creates, in a Dual context, the Backward Star and the + // Forward Star of every node + + for( nodeDType *node = nodesD ; node != stopNodesD ; node++) { + // initialize nodes + node->firstBs = NULL; + node->firstFs = NULL; + node->numArcs = 0; + } + + dummyRootD->firstBs = NULL; + dummyRootD->firstFs = NULL; + dummyRootD->numArcs = 0; + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) { + // initialize real arcs + arc->nextFs = ( arc->tail )->firstFs; + ( arc->tail )->firstFs = arc; + arc->nextBs = ( arc->head )->firstBs; + ( arc->head )->firstBs = arc; + ( arc->tail )->numArcs++; + ( arc->head )->numArcs++; + } + + ResetWhenInT2(); + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::PrimalSimplex( void ) +{ + #if( UNIPI_PRIMAL_INITIAL_SHOW == 0 ) + #if( UNIPI_PRIMAL_ITER_SHOW == 0 ) + #if( UNIPI_PRIMAL_FINAL_SHOW == 0 ) + //cout << endl; + #endif + #endif + #endif + #if( UNIPI_PRIMAL_INITIAL_SHOW ) + cout << endl; + for( int t = 0; t < 3; t++ ) + cout << "\t"; + cout << "PRIMALE MCFSimplex: ARCHI E NODI ALL' INIZIO" << endl; + ShowSituation( 3 ); + #endif + #if( QUADRATICCOST ) + #if( LIMITATEPRECISION ) + foValue = GetFO(); + int cont = 0; + #endif + #endif + + status = kUnSolved; + if( pricingRule != kCandidateListPivot ) + arcToStartP = arcsP; + + iterator = 0; // setting the initial arc for the Dantzig or First Elibigle Rule + + arcPType *enteringArc; + arcPType *leavingArc; + if( pricingRule == kCandidateListPivot ) + InitializePrimalCandidateList(); + + while( status == kUnSolved ) { + iterator++; + switch( pricingRule ) { + case( kDantzig ): + enteringArc = RuleDantzig(); + break; + case( kFirstEligibleArc ): + enteringArc = PRuleFirstEligibleArc(); + break; + case( kCandidateListPivot ): + enteringArc = RulePrimalCandidateListPivot(); + break; + } + + #if( QUADRATICCOST ) + #if( LIMITATEPRECISION ) + /* In the quadratic case with LIMITATEPRECISION == 1, the entering arcs + are selected according to a thresold. + This thresold is definited according to the old f.o. value. + If Primal Simplex doesn't find an entering arc, it calculates again + the f.o. value, and try again. */ + if( enteringArc == NULL ) { + foValue = GetFO(); + switch( pricingRule ) { + case( kDantzig ): + enteringArc = RuleDantzig(); + break; + case( kFirstEligibleArc ): + enteringArc = PRuleFirstEligibleArc(); + break; + case( kCandidateListPivot ): + enteringArc = RulePrimalCandidateListPivot(); + break; + } + } + #endif + #endif + + if( pricingRule != kCandidateListPivot ) { + // in every iteration the algorithm changes the initial arc for + // Dantzig and First Eligible Rule. + arcToStartP++; + if( arcToStartP == stopArcsP ) + arcToStartP = arcsP; + } + + if( enteringArc ) { + arcPType *arc; + nodePType *k1; + nodePType *k2; + /* If the reduced cost of the entering arc is > 0, the Primal Simplex + pushes flow in the cycle determinated by T and entering arc for decreases + flow in the entering arc: in the linear case entering arc's flow goes to 0, + in the quadratic case it decreases while it's possibile. + If the reduced cost of the entering arc is < 0, the Primal Simplex pushes + flow in the cycle determinated by T and entering arc for increases flow + in the entering arc: in the linear case entering arc's flow goes to upper + bound, in the quadratic case it increases while it's possibile. */ + + #if( QUADRATICCOST ) + FONumber t; + FONumber theta; + FONumber deltaFO; + FNumber theta2; + CNumber Q = ( enteringArc->tail )->sumQuadratic + + ( enteringArc->head )->sumQuadratic + enteringArc->quadraticCost; + // Q is the sum of the quadratic coefficient in the cycle determinated by T + // and entering arc. + FONumber rc = ReductCost( enteringArc ); + if( ETZ( Q, EpsCst ) ) + theta = Inf<FNumber>(); // This value will be certainly decreased + else + theta = ABS( rc / Q ); + // This is the best theta value (with best f.o. value decrease) + + leavingArc = enteringArc; + nodePType *cycleRoot; // The radix of the cycle determinated by T and entering arc. + if( GTZ( rc , EpsCst ) ) { + #else + FNumber t; + FNumber theta; + if( enteringArc->ident == AT_UPPER ) { + #endif + /* Primal Simplex increases or decreases entering arc's flow. + "theta" is a positive value. + For this reason the algorithm uses two nodes ("k1" and "k2") to push + flow ("theta") from "k1" to "k2". + According to entering arc's reduct cost, the algorithm determinates + "k1" and "k2" */ + + k1 = enteringArc->head; + k2 = enteringArc->tail; + #if( QUADRATICCOST ) + theta = min( theta , enteringArc->flow ); + // The best value for theta is compared with the entering arc's bound + theta2 = - theta; + #else + theta = enteringArc->flow; + #endif + } + else { + k1 = enteringArc->tail; + k2 = enteringArc->head; + #if( QUADRATICCOST ) + theta = min( theta , enteringArc->upper - enteringArc->flow ); + // The best value for theta is compared with the entering arc's bound + theta2 = theta; + #else + theta = enteringArc->upper - enteringArc->flow; + #endif + } + + nodePType *memK1 = k1; + nodePType *memK2 = k2; + leavingArc = NULL; + #if( QUADRATICCOST ) + #if( LIMITATEPRECISION ) + deltaFO = rc * theta2 + Q * theta2 * theta2 / 2; + #endif + bool leavingReducesFlow = GTZ( rc , EpsCst ); + #else + bool leavingReducesFlow = GTZ( ReductCost( enteringArc ) , EpsCst ); + #endif + // Compute "theta", find outgoing arc and "root" of the cycle + bool leave; + // Actual "theta" is compared with the bounds of the other cycle's arcs + while( k1 != k2 ) { + if( k1->subTreeLevel > k2->subTreeLevel ) { + arc = k1->enteringTArc; + if( arc->tail != k1 ) { + t = arc->upper - arc->flow; + leave = false; + } + else { + t = arc->flow; + leave = true; + } + + if( t < theta ) { + // The algorithm has found a possible leaving arc + theta = t; + leavingArc = arc; + leavingReducesFlow = leave; + // If "leavingReducesFlow" == true, if this arc is selected to exit the base, + // it decreases its flow + } + + k1 = Father( k1 , arc ); + } + else { + arc = k2->enteringTArc; + if( arc->tail == k2 ) { + t = arc->upper - arc->flow; + leave = false; + } + else { + t = arc->flow; + leave = true; + } + + if( t <= theta ) { + // The algorithm has found a possible leaving arc + theta = t; + leavingArc = arc; + leavingReducesFlow = leave; + // If "leavingReducesFlow" == true, if this arc is selected to exit the base, + // it decreases its flow + } + + k2 = Father(k2, arc); + } + } + + #if( QUADRATICCOST ) + cycleRoot = k1; + #endif + + if( leavingArc == NULL ) + leavingArc = enteringArc; + + if( theta >= Inf<FNumber>() ) { + status = kUnbounded; + break; + } + + // Update flow with "theta" + k1 = memK1; + k2 = memK2; + #if( QUADRATICCOST ) + if( enteringArc->tail == k1 ) + theta2 = theta; + else + theta2 = -theta; + + // "theta" is a positive value in every case. + // "theta2" is the real theta value according to the real + // direction of the entering arc + #if( LIMITATEPRECISION ) + deltaFO = rc * theta2 + Q * theta2 * theta2 / 2; + // The decrease of the f.o. value in the quadratic case + #endif + #endif + + if( ! ETZ(theta , EpsFlw ) ) { + if( enteringArc->tail == k1 ) + enteringArc->flow = enteringArc->flow + theta; + else + enteringArc->flow = enteringArc->flow - theta; + + while( k1 != k2 ) { + if( k1->subTreeLevel > k2->subTreeLevel ) { + arc = k1->enteringTArc; + if( arc->tail != k1 ) + arc->flow = arc->flow + theta; + else + arc->flow = arc->flow - theta; + + k1 = Father(k1, k1->enteringTArc); + } + else { + arc = k2->enteringTArc; + if( arc->tail == k2 ) + arc->flow = arc->flow + theta; + else + arc->flow = arc->flow - theta; + + k2 = Father(k2, k2->enteringTArc); + } + } + } + + if( enteringArc != leavingArc ) { + bool leavingBringFlowInT2 = ( leavingReducesFlow == + ( ( leavingArc->tail )->subTreeLevel > ( leavingArc->head )->subTreeLevel ) ); + // "leavingBringFlowInT2" == true if leaving arc brings flow to the subtree T2 + if( leavingBringFlowInT2 == ( memK1 == enteringArc->tail ) ) { + k2 = enteringArc->tail; + k1 = enteringArc->head; + } + else { + k2 = enteringArc->head; + k1 = enteringArc->tail; + } + } + + #if( QUADRATICCOST == 0 ) + if( leavingReducesFlow ) + leavingArc->ident = AT_LOWER; + else + leavingArc->ident = AT_UPPER; + + if( leavingArc != enteringArc ) { + enteringArc->ident = BASIC; + nodePType *h1; + nodePType *h2; + // "h1" is the node in the leaving arc with smaller tree's level + if( ( leavingArc->tail )->subTreeLevel < ( leavingArc->head )->subTreeLevel ) { + h1 = leavingArc->tail; + h2 = leavingArc->head; + } + else { + h1 = leavingArc->head; + h2 = leavingArc->tail; + } + + UpdateT(leavingArc, enteringArc, h1, h2, k1, k2); + // Update potential of the subtree T2 + k2 = enteringArc->head; + CNumber delta = ReductCost(enteringArc); + if( ( enteringArc->tail )->subTreeLevel > ( enteringArc->head )->subTreeLevel ) { + delta = -delta; + k2 = enteringArc->tail; + } + + AddPotential( k2 , delta ); + // In the linear case Primal Simplex only updates the potential of the nodes of + // subtree T2 + } + #else + if( leavingArc != enteringArc ) { + nodePType *h1; + nodePType *h2; + // "h1" is the node in the leaving arc with smaller tree's level + if( ( leavingArc->tail )->subTreeLevel < + ( leavingArc->head )->subTreeLevel ) { + h1 = leavingArc->tail; + h2 = leavingArc->head; + } + else { + h1 = leavingArc->head; + h2 = leavingArc->tail; + } + + // Update the basic tree + UpdateT( leavingArc , enteringArc , h1 , h2 , k1 , k2 ); + } + + #if( OPTQUADRATIC ) + nodePType *h1; + nodePType *h2; + if( ( leavingArc->tail )->subTreeLevel < + ( leavingArc->head )->subTreeLevel ) { + h1 = leavingArc->tail; + h2 = leavingArc->head; + } + else { + h1 = leavingArc->head; + h2 = leavingArc->tail; + } + + nodePType *node = h1; + nodePType *updateNode = h1; + if( h1 == cycleRoot ) + ComputePotential( cycleRoot ); + else { + while( node != cycleRoot ) { + arcPType *entArc = node->enteringTArc; + if( ! ETZ( entArc->quadraticCost , EpsCst ) ) + updateNode = node; + + node = Father( node , entArc ); + } + + ComputePotential( updateNode ); + node = h2; + updateNode = h2; + while( node != cycleRoot ) { + arcPType *entArc = node->enteringTArc; + if( ! ETZ( entArc->quadraticCost , EpsCst ) ) + updateNode = node; + + node = Father( node , entArc ); + } + + ComputePotential( updateNode ); + } + #else + // Update the potential of the node "manually" + ComputePotential( cycleRoot ); + #endif + + #if( LIMITATEPRECISION ) + cont = cont + 1; + if( cont == recomputeFOLimits ) { + cont = 0; + foValue = GetFO(); // Calculate f.o. value manually + } + else + foValue = foValue + deltaFO; + // Calculate the f.o. value with the estimated decrease + #endif + #endif + } + else { + status = kOK; + // If one of dummy arcs has flow bigger than 0, the solution is unfeasible. + for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) + if( GTZ( arc->flow , EpsFlw ) ) + status = kUnfeasible; + } + + if( ( status == kUnSolved ) && MaxTime && MCFt ) { + double tu, ts; + TimeMCF( tu , ts ); + if( MaxTime < tu + ts ) + status = kStopped; + } + + if( ( status == kUnSolved ) && MaxIter) + if( MaxIter < (int) iterator ) + status = kStopped; + + #if( UNIPI_PRIMAL_ITER_SHOW ) + int it = (int) iterator; + if( it % UNIPI_PRIMAL_ITER_SHOW == 0 ) { + cout << endl; + for( int t = 0; t < 3; t++ ) + cout << "\t"; + cout << "PRIMALE MCFSimplex: ARCHI E NODI ALLA " << it << " ITERAZIONE" << endl; + ShowSituation( 3 ); + #if( FOSHOW ) + if( (int) iterator % FOSHOW == 0 ) + clog << "Iteration = " << iterator << " of = " + #if( LIMITATEPRECISION && QUADRATICCOST ) + << foValue + #else + << GetFO() + #endif + << endl; + #endif + } + #endif + } + + #if( UNIPI_PRIMAL_FINAL_SHOW ) + cout << endl; + for( int t = 0; t < 3; t++ ) + cout << "\t"; + cout << "PRIMALE UniPi: ARCHI E NODI ALLA FINE" << endl; + ShowSituation( 3 ); + #endif + + } // end( PrimalSimplex ) + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::DualSimplex( void ) +{ + #if( UNIPI_PRIMAL_INITIAL_SHOW == 0 ) + #if( UNIPI_PRIMAL_ITER_SHOW == 0 ) + #if( UNIPI_PRIMAL_FINAL_SHOW == 0 ) + cout << endl; + #endif + #endif + #endif + #if( UNIPI_DUAL_INITIAL_SHOW ) + cout << endl; + for( int t = 0; t < 3; t++ ) + cout << "\t"; + cout << "DUALE MCFSimplex: ARCHI E NODI ALL' INIZIO" << endl; + ShowSituation( 3 ); + #endif + + if( pricingRule != kCandidateListPivot ) + arcToStartD = arcsD; + + iterator = 0; + arcDType *enteringArc; + arcDType *leavingArc; + if( pricingRule == kCandidateListPivot ) + InitializeDualCandidateList(); + + status = kUnSolved; + while( status == kUnSolved ) { + iterator++; + if( iterator == Inf<iteratorType>() ) { + ResetWhenInT2(); // Restore to 0 every nodes' field "whenInT2" + iterator = 1; + } + + switch( pricingRule ) { + case( kDantzig ): + leavingArc = DRuleFirstEligibleArc(); + break; // si esegue cmq FEA + case( kFirstEligibleArc ): + leavingArc = DRuleFirstEligibleArc(); + break; + case( kCandidateListPivot ): + leavingArc = RuleDualCandidateListPivot(); + break; + } + + if( pricingRule != kCandidateListPivot ) { + arcToStartD++; + if( arcToStartD == stopArcsD ) + arcToStartD = dummyArcsD; + + if( arcToStartD == stopDummyD ) + arcToStartD = arcsD; + // Setting the initial arc for the Dantzig or First Elibigle Rule + } + + if( leavingArc ) { + bool leavingArcInL = false; + bool leavingArcFromT1toT2; + if( LTZ( leavingArc->flow , EpsFlw ) ) + leavingArcInL = true; + + nodeDType *h1; + nodeDType *h2; + if( ( leavingArc->tail )->subTreeLevel < + ( leavingArc->head )->subTreeLevel ) { + h1 = leavingArc->tail; + h2 = leavingArc->head; + leavingArcFromT1toT2 = true; + } + else { + h1 = leavingArc->head; + h2 = leavingArc->tail; + leavingArcFromT1toT2 = false; + } + + Index numOfT2Arcs = 0; + int level = h2->subTreeLevel; + nodeDType *node = h2; + node->whenInT2 = iterator; + nodeDType *lastNodeOfT2 = h2; + numOfT2Arcs = node->numArcs; + while( node->nextInT && ( ( node->nextInT )->subTreeLevel > level ) ) { + node = node->nextInT; + lastNodeOfT2 = node; + numOfT2Arcs = numOfT2Arcs + node->numArcs; + node->whenInT2 = iterator; + } + + /* The Dual Simplex has determinated the leaving arc, and so the + subtrees T1 and T2. Dual Simplex scans T2 to fix the fields "whenInT2" + of T2's nodes to the iteration value, and counts the entering and + outgoing arcs from these nodes. According to this number, it decides + to scan the Backward and Forward of the subtree (T1 or T2) with the + minor number of entering/outgoing arcs from its nodes. */ + enteringArc = NULL; + bool lv = ( leavingArcFromT1toT2 == leavingArcInL ); + CNumber maxRc = Inf<CNumber>(); + //Search arc in the Forward Star and Backward Star of nodes of T1 + if( numOfT2Arcs > m ) { + // Dual Simplex starts from the node which follows the dummy root. + node = dummyRootD->nextInT; + bool fine = false; + while( fine == false ) { + /* If node is the root of subtree T2, Dual Simplex jumps to the node + (if exists) which follows the last node of T2 */ + if( node == h2 ) + if( lastNodeOfT2->nextInT ) + node = lastNodeOfT2->nextInT; + else + break; + + // Search arc in the Backward Star of nodes of T1 + arcDType *arc = node->firstBs; + while( arc ) { + if( ( arc->tail )->whenInT2 == iterator ) { + // Evaluate each arc from T2 to T1 which isn't in T + if( arc->ident == AT_LOWER ) { + if( lv ) { + CNumber rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst) ) { + fine = true; + break; + } + } + } + } + + if( arc->ident == AT_UPPER ) { + if( ! lv ) { + CNumber rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + + if( ETZ( rc , EpsCst ) ) { + fine = true; + break; + } + } + } + } + } + + arc = arc->nextBs; + } + + // Search arc in the Forward Star of nodes of T1 + arc = node->firstFs; + while( arc ) { + if( ( arc->head )->whenInT2 == iterator ) { + // Evaluate each arc from T1 to T2 which isn't in T + if( arc->ident == AT_LOWER ) { + if( ! lv ) { + CNumber rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst ) ) { + fine = true; + break; + } + } + } + } + + if( arc->ident == AT_UPPER ) { + if( lv ) { + CNumber rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst ) ) { + fine = true; + break; + } + } + } + } + } + + arc = arc->nextFs; + } + + node = node->nextInT; + if( node == NULL ) + fine = true; + } + } + // Search arc in the Forward Star and Backward Star of nodes of T2 + else { + node = h2; + bool fine = false; + while( fine == false ) { + // Search arc in the Backward Star of nodes of T2 + arcDType *arc = node->firstBs; + CNumber rc; + while( arc ) { + if( ( arc->tail )->whenInT2 != iterator ) { + // Evaluate each arc from T1 to T2 which isn't in T + if( arc->ident == AT_LOWER ) { + if( ! lv ) { + rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst ) ) { + fine = true; + break; + } + } + } + } + + if( arc->ident == AT_UPPER ) { + if( lv ) { + rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst ) ) { + fine = true; + break; + } + } + } + } + } + + arc = arc->nextBs; + } + + // Search arc in the Forward Star of nodes of T2 + arc = node->firstFs; + while( arc ) { + if( ( arc->head )->whenInT2 != iterator ) { + // Evaluate each arc from T2 to T1 which isn't in T + if( arc->ident == AT_LOWER ) { + if( lv ) { + rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst ) ) { + fine = true; + break; + } + } + } + } + + if( arc->ident == AT_UPPER ) { + if( ! lv ) { + rc = ABS( ReductCost( arc ) ); + if( LT( rc , maxRc , EpsCst ) ) { + enteringArc = arc; + maxRc = rc; + /* If arc is appropriate to enter in T and its reduct cost is 0, + search is ended: this is the arc which enters in T */ + if( ETZ( rc , EpsCst) ) { + fine = true; + break; + } + } + } + } + } + + arc = arc->nextFs; + } + + if( node == lastNodeOfT2 ) + fine = true; + else + node = node->nextInT; + + } + } + + if( enteringArc ) { + FNumber theta = -leavingArc->flow; + if( GTZ( leavingArc->flow , EpsFlw ) ) + theta = leavingArc->flow - leavingArc->upper; + // Initial value of theta is the infeasibility of the leaving arc + + FNumber t; + nodeDType *k1; + nodeDType *k2; + /* if entering arc is in U, Dual Simplex pushs flow in the cycle + determinated by T and entering arc for decreases flow in the entering arc: + if entering arc is in L, Dual Simplex pushs flow in the cycle + determinated by T and entering arc for increases flow in the entering arc: + Dual Simplex increases or decreases entering arc's flow. + theta is a positive value. + For this reason the algorithm uses two nodes (k1 and k2) to push flow + (theta) from k1 to k2. According to entering arc's reduct cost, the + algorithm determinates k1 and k2 */ + + if( enteringArc->ident == AT_UPPER ) { + k1 = enteringArc->head; + k2 = enteringArc->tail; + } + else { + k1 = enteringArc->tail; + k2 = enteringArc->head; + } + + nodeDType *memK1 = k1; + nodeDType *memK2 = k2; + arcDType *arc; + k1 = memK1; + k2 = memK2; + // Update the flow + while( k1 != k2 ) { + if( k1->subTreeLevel > k2->subTreeLevel ) { + arc = k1->enteringTArc; + if( arc->tail != k1 ) + arc->flow = arc->flow + theta; + else + arc->flow = arc->flow - theta; + + k1 = Father(k1, k1->enteringTArc); + } + else { + arc = k2->enteringTArc; + if( arc->tail == k2 ) + arc->flow = arc->flow + theta; + else + arc->flow = arc->flow - theta; + + k2 = Father( k2 , k2->enteringTArc ); + } + } + + if(leavingArcInL ) + leavingArc->ident = AT_LOWER; + else + leavingArc->ident = AT_UPPER; + + bool leavingBringFlowInT2 = ( leavingArcInL == + ( ( leavingArc->tail )->subTreeLevel > + ( leavingArc->head )->subTreeLevel ) ); + // leavingBringFlowInT2 == true if leaving arc brings flow to the subtree T2 + if( leavingBringFlowInT2 != ( memK1 == enteringArc->tail ) ) { + k2 = enteringArc->tail; + k1 = enteringArc->head; + } + else { + k2 = enteringArc->head; + k1 = enteringArc->tail; + } + + if( enteringArc->ident == AT_LOWER ) + enteringArc->flow = enteringArc->flow + theta; + else + enteringArc->flow = enteringArc->flow - theta; + + enteringArc->ident = BASIC; + UpdateT( leavingArc , enteringArc , h1 , h2 , k1 , k2 ); + // update potential of the subtree T2 + k2 = enteringArc->head; + CNumber delta = ReductCost( enteringArc ); + if( ( enteringArc->tail) ->subTreeLevel > + ( enteringArc->head )->subTreeLevel ) { + delta = -delta; + k2 = enteringArc->tail; + } + + // Dual Simplex only updates the potential of the T2's nodes + AddPotential( k2 , delta ); + } + else + status = kUnfeasible; + /* If Dual Simplex finds a leaving arc but it doesn't find an entering arc, + the algorithm stops. At this point Dual Simplex has an unfeasible primal + solution. */ + } + else { + status = kOK; + // If one of dummy arcs has flow different than 0, the solution is unfeasible. + for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ ) + if( ! ETZ( arc->flow , EpsFlw ) ) { + status = kUnfeasible; + break; + } + } + + if( ( status == kUnSolved ) && MaxTime ) { + double tu, ts; + TimeMCF( tu , ts ); + if( MaxTime < tu + ts ) + status = kStopped; + } + + if( ( status == kUnSolved ) && MaxIter && MCFt ) + if( MaxIter < (int) iterator ) + status = kStopped; + + #if( UNIPI_DUAL_ITER_SHOW ) + if( (int) iterator % UNIPI_DUAL_ITER_SHOW == 0 ) { + cout << endl; + for( int t = 0; t < 3; t++ ) + cout << "\t"; + cout << "DUALE MCFSimplex: ARCHI E NODI ALLA " << iterator << " ITERAZIONE" + << endl; + ShowSituation( 3 ); + #if( FOSHOW ) + cout << "of = " << GetFO() << endl; + #endif + } + #endif + } + + #if( UNIPI_DUAL_ITER_SHOW ) + int it = (int) iterator; + if( it % UNIPI_DUAL_ITER_SHOW == 0 ) { + cout << endl; + for( int t = 0; t < 3; t++ ) + cout << "\t"; + cout << "DUALE MCFSimplex: ARCHI E NODI ALLA " << iterator << " ITERAZIONE" + << endl; + Showsituation( 3 ); + } + #endif + + } // end( DualSimplex ) + +/*--------------------------------------------------------------------------*/ + +template<class N, class A> +void MCFSimplex::UpdateT( A *h , A *k , N *h1 , N *h2 , N *k1 , N *k2 ) +{ + /* In subtree T2 there is a path from node h2 (deepest node of the leaving + arc h and root of T2) to node k2 (deepest node of the leaving arc h and + coming root of T2). With the update of T, this path will be overturned: + node k2 will become the root of T2... + The subtree T2 must be reordered and the field "subTreeLevel", which + represents the depth in T of every node, of every T2's nodes is changed. + Variable delta represents the increase of "subTreeLevel" value for node + k2 and its descendants: probably this value is a negative value. */ + + int delta = (k1->subTreeLevel) + 1 - (k2->subTreeLevel); + N *root = k2; + N *dad; + + /*To reorder T2, the method analyses every nodes of the path h2->k2, starting + from k2. For every node, it moves the node's subtree from its original + position to a new appropriate position. In particular k2 and its subtree + (k2 is the new root of T2, so the first nodes of the new T2) will be moved + next to node k1 (new father of k2), the next subtree will be moved beside + the last node of k2's subtree.... + "previousNode" represents the node where the new subtree will be moved + beside in this iterative action. At the start "previousNode" is the node + k1 (T2 will be at the right of k1). */ + + N *previousNode = k1; + N *lastNode; + /* "arc1" is the entering arc in T (passed by parameters). + For every analysed node of path h2->k2, the method changes + "enteringTArc" but it must remember the old arc, which will be the + "enteringTArc" of the next analysed node. + At the start "arc1" is k (the new "enteringTArc" of k2). */ + + A *arc1 = k; + A *arc2; + bool fine = false; + while( fine == false ) { + // If analysed node is h2, this is the last iteration + if( root == h2 ) + fine = true; + + dad = Father( root , root->enteringTArc ); + // Cut the root's subtree from T and update the "subLevelTree" of its nodes + lastNode = CutAndUpdateSubtree( root , delta ); + // Paste the root's subtree in the right position; + PasteSubtree( root , lastNode , previousNode ); + // In the next iteration the subtree will be moved beside the last node of + // the actual analysed subtree. + previousNode = lastNode; + /* The increase of the subtree in the next iteration is different from + the actual increase: in particual the increase increases itself (+2 + at every iteration). */ + delta = delta + 2; + /* A this point "enteringTArc" of actual root is stored in "arc2" and + changed; then "arc1" and "root" are changed. */ + arc2 = root->enteringTArc; + root->enteringTArc = arc1; + arc1 = arc2; + root = dad; + } + } + +/*--------------------------------------------------------------------------*/ + +template<class N> +N* MCFSimplex::CutAndUpdateSubtree( N *root , int delta ) +{ + int level = root->subTreeLevel; + N *node = root; + // The root of this subtree is passed by parameters, the last node is searched. + while ( ( node->nextInT ) && ( ( node->nextInT )->subTreeLevel > level ) ) { + node = node->nextInT; + // The "subTreeLevel" of every nodes of subtree is updated + node->subTreeLevel = node->subTreeLevel + delta; + } + + root->subTreeLevel = root->subTreeLevel + delta; + /* The 2 neighbouring nodes of the subtree (the node at the left of the root + and the node at the right of the last node) is joined. */ + + if( root->prevInT ) + ( root->prevInT )->nextInT = node->nextInT; + if( node->nextInT ) + ( node->nextInT )->prevInT = root->prevInT; + + return( node ); // the method returns the last node of the subtree + } + +/*--------------------------------------------------------------------------*/ + +template<class N> +void MCFSimplex::PasteSubtree( N *root , N *lastNode , N *previousNode ) +{ + /* The method inserts subtree ("root" and "lastNode" are the extremity of the + subtree) after "previousNode". The method starts to identify the next node + of "previousNode" ("nextNode"), so it joins "root" with "previousNode" and + "lastNode" with "nextNode" (if exists). */ + + N *nextNode = previousNode->nextInT; + root->prevInT = previousNode; + previousNode->nextInT = root; + lastNode->nextInT = nextNode; + if( nextNode ) + nextNode->prevInT = lastNode; + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcPType* MCFSimplex::RuleDantzig( void ) +{ + arcPType *arc = arcToStartP; + arcPType *enteringArc = NULL; + #if( QUADRATICCOST ) + /* In the quadratic case used type for reduct cost is FONumber. + Value "lim" is the fixed thresold for the decrease of the f.o. value */ + FONumber lim = EpsOpt * foValue / n; + FONumber RC; + FONumber maxValue = 0; + #else + CNumber RC; + CNumber maxValue = 0; + #endif + + do { + // The method analyses every arc + #if( QUADRATICCOST ) + RC = ReductCost( arc ); + FNumber theta; + /* If reduct cost of arc is lower than 0, the flow of the arc must increase. + If reduct cost of arc is bigger than 0, the flow of the arc must decrease. + "theta" is the difference between lower (upper) bound and the actual flow. + */ + + if( LTZ( RC , EpsCst ) ) + theta = arc->upper - arc->flow; + + if( GTZ( RC , EpsCst ) ) + theta = -arc->flow; + + // If it's possible to increase (or decrease) the flow in this arc + if( ! ETZ( theta , EpsFlw ) ) { + /* "Q" is the sum of the quadratic coefficient of the arc belonging the + T path from tail's arc to head's arc + "Q" is always bigger than 0 or equals to 0. + If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow + with the best decrease of f.o. value. + - RC/ Q must be compare with "theta" to avoid that the best increase + (decrease) of the flow violates the bounds of the arc. + This confront determines "theta". */ + + CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic + + arc->quadraticCost; + + if( GTZ( Q , EpsCst ) ) + if( GTZ( theta , EpsFlw ) ) + theta = min( theta , - RC / Q ); + else + theta = max( theta , - RC / Q ); + + /* Calculate the estimate decrease of the f.o. value if this arc is + selected and flow is increased (decreased) by "theta" */ + + CNumber deltaFO = RC * theta + Q * theta * theta / 2; + // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than + // old decrease value, arc is the best arc. + + if( deltaFO < maxValue ) { + maxValue = deltaFO; + enteringArc = arc; + } + } + #else + if( arc->ident > BASIC ) { + RC = ReductCost( arc ); + + if( ( LTZ( RC , EpsCst ) && ( arc->ident == AT_LOWER ) ) || + ( GTZ( RC , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) { + + if( RC < 0 ) + RC = -RC; + + if( RC > maxValue ) { + maxValue = RC; + enteringArc = arc; + } + } + } + #endif + + arc++; + if( arc == stopArcsP ) + arc = arcsP; + + } while( arc != arcToStartP ); + + #if( ( LIMITATEPRECISION ) && ( QUADRATICCOST ) ) + if( -maxValue <= lim ) + enteringArc = NULL; + #endif + + return( enteringArc ); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcPType* MCFSimplex::PRuleFirstEligibleArc( void ) +{ + arcPType *arc = arcToStartP; + arcPType *enteringArc = NULL; + #if( QUADRATICCOST ) + FONumber RC; + #else + CNumber RC; + #endif + + do { + #if( QUADRATICCOST ) + // In this method the "decrease f.o. value" criterion is not used. + RC = ReductCost( arc ); + if( ( LTZ( RC , EpsCst ) && LT( arc->flow , arc->upper , EpsFlw ) ) || + ( GTZ( RC , EpsCst ) && GTZ( arc->flow , EpsFlw ) ) ) + enteringArc = arc; + #else + if( arc->ident > BASIC ) { + RC = ReductCost( arc ); + if( ( LTZ( RC , EpsCst ) && ( arc->ident == AT_LOWER ) ) || + ( GTZ( RC , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) + enteringArc = arc; + } + #endif + + arc++; + if( arc == stopArcsP ) + arc = dummyArcsP; + if( arc == stopDummyP ) + arc = arcsP; + + } while( ( enteringArc == NULL ) && ( arc != arcToStartP ) ); + + return( enteringArc ); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcDType* MCFSimplex::DRuleFirstEligibleArc( void ) +{ + arcDType *arc = arcToStartD; + arcDType *leavingArc = NULL; + do { + if( GT( arc->flow , arc->upper , EpsFlw ) || LTZ( arc->flow , EpsFlw ) ) + leavingArc = arc; + + arc++; + if( arc == stopArcsD ) + arc = dummyArcsD; + if( arc == stopDummyD ) + arc = arcsD; + + } while( ( leavingArc == NULL ) && ( arc != arcToStartD ) ); + + return(leavingArc); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcPType* MCFSimplex::RulePrimalCandidateListPivot( void ) +{ + Index next = 0; + Index i; + Index minimeValue; + if( hotListSize < tempCandidateListSize ) + minimeValue = hotListSize; + else + minimeValue = tempCandidateListSize; + + #if( QUADRATICCOST ) + // Check if the left arcs in the list continue to violate the dual condition + for( i = 2 ; i <= minimeValue ; i++ ) { + arcPType *arc = candP[ i ].arc; + FONumber red_cost = ReductCost( arc ); + FNumber theta = 0; + /* If reduct cost of arc is lower than 0, the flow of the arc must increase. + If reduct cost of arc is bigger than 0, the flow of the arc must decrease. + "theta" is the difference between lower (upper) bound and the actual flow. + */ + + if( LTZ( red_cost , EpsCst ) ) + theta = arc->upper - arc->flow; + else + if( GTZ( red_cost , EpsCst ) ) + theta = - arc->flow; + + // if it's possible to increase (or decrease) the flow in this arc + if( theta != 0 ) { + /* "Q" is the sum of the quadratic coefficient of the arc belonging the T + path from tail's arc to head's arc + "Q" is always bigger than 0 or equals to 0. + If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow + with the best decrease of f.o. value. + - RC/ Q must be compare with "theta" to avoid that the best increase + (decrease) of the flow violates the bounds of the arc. + This confront determines "theta". */ + + CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic + + arc->quadraticCost; + + if( GTZ( Q , EpsCst ) ) + if( GTZ( theta , EpsFlw ) ) + theta = min( theta , - red_cost / Q ); + else + theta = max( theta , - red_cost / Q ); + + /* Calculate the estimate decrease of the f.o. value if this arc is + selected and flow is increased (decreased) by "theta" */ + + CNumber deltaFO = red_cost * theta + Q * theta * theta / 2; + #if( LIMITATEPRECISION ) + // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than + // old decrease value, arc is the best arc. + if( - deltaFO > ( EpsOpt * foValue / n ) ) { + #else + if( LTZ( deltaFO , EpsCst ) ) { + #endif + next++; + candP[ next ].arc = arc; + candP[ next ].absRC = -deltaFO; + } + } + } + + tempCandidateListSize = next; + Index oldGroupPos = groupPos; + // Search other arcs to fill the list + do { + arcPType *arc; + for( arc = arcsP + groupPos ; arc < stopArcsP ; arc += numGroup ) { + FONumber red_cost = ReductCost( arc ); + FNumber theta = 0; + /* If reduced cost is lower than 0, the flow of the arc must increase. + If reduced cost is larger than 0, the flow of the arc must decrease. + "theta" is the difference between lower (upper) bound and the actual + flow. */ + + if( LTZ( red_cost , EpsCst ) ) + theta = arc->upper - arc->flow; + else + if( GTZ( red_cost , EpsCst ) ) + theta = - arc->flow; + + // if it's possible to increase (or decrease) the flow in this arc + if( theta != 0 ) { + /* "Q" is the sum of the quadratic coefficient of the arc belonging the T + path from tail's arc to head's arc + "Q" is always bigger than 0 or equals to 0. + If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow + with the best decrease of f.o. value. + - RC/ Q must be compare with "theta" to avoid that the best increase + (decrease) of the flow violates the bounds of the arc. + This confront determines "theta". */ + CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic + + arc->quadraticCost; + + if( GTZ( Q , EpsCst ) ) + if( GTZ( theta , EpsFlw ) ) + theta = min( theta , - red_cost / Q ); + else + theta = max( theta , - red_cost / Q ); + + /* Calculate the estimate decrease of the f.o. value if this arc is + selected and flow is increased (decreased) by "theta" */ + + CNumber deltaFO = red_cost * theta + Q * theta * theta / 2; + #if( LIMITATEPRECISION ) + // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than + // old decrease value, arc is the best arc. + if( -deltaFO > ( EpsOpt * foValue / n ) ) { + #else + if( LTZ( deltaFO , EpsCst ) ) { + #endif + tempCandidateListSize++; + candP[ tempCandidateListSize ].arc = arc; + candP[ tempCandidateListSize ].absRC = -deltaFO; + } + } + } + + groupPos++; + if( groupPos == numGroup ) + groupPos = 0; + + } while( ( tempCandidateListSize < hotListSize ) && + ( groupPos != oldGroupPos ) ); + #else + // Check if the left arcs in the list continue to violate the dual condition + for( i = 2 ; i <= minimeValue ; i++ ) { + arcPType *arc = candP[i].arc; + CNumber red_cost = ReductCost( arc ); + + if( ( LTZ( red_cost , EpsCst ) && ( arc->ident == AT_LOWER ) ) || + ( GTZ( red_cost , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) { + next++; + candP[ next ].arc = arc; + candP[ next ].absRC = ABS( red_cost ); + } + } + + tempCandidateListSize = next; + Index oldGroupPos = groupPos; + // Search other arcs to fill the list + do { + arcPType *arc; + for( arc = arcsP + groupPos ; arc < stopArcsP ; arc += numGroup ) { + if( arc->ident == AT_LOWER ) { + CNumber red_cost = ReductCost( arc ); + if( LTZ( red_cost , EpsCst ) ) { + tempCandidateListSize++; + candP[ tempCandidateListSize ].arc = arc; + candP[ tempCandidateListSize ].absRC = ABS( red_cost ); + } + } + else + if( arc->ident == AT_UPPER ) { + CNumber red_cost = ReductCost( arc ); + if( GTZ( red_cost , EpsCst ) ) { + tempCandidateListSize++; + candP[ tempCandidateListSize ].arc = arc; + candP[ tempCandidateListSize ].absRC = ABS( red_cost ); + } + } + } + + groupPos++; + if( groupPos == numGroup ) + groupPos = 0; + + } while( ( tempCandidateListSize < hotListSize ) && ( groupPos != oldGroupPos ) ); + #endif + + if( tempCandidateListSize ) { + SortPrimalCandidateList( 1 , tempCandidateListSize ); + return( candP[ 1 ].arc ); + } + else + return( NULL ); + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFSimplex::InitializePrimalCandidateList( void ) +{ + numGroup = ( ( m - 1 ) / numCandidateList ) + 1; + groupPos = 0; + tempCandidateListSize = 0; + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFSimplex::SortPrimalCandidateList( Index min , Index max ) +{ + Index left = min; + Index right = max; + #if( QUADRATICCOST ) + FONumber cut = candP[ ( left + right ) / 2 ].absRC; + #else + CNumber cut = candP[ ( left + right ) / 2 ].absRC; + #endif + do { + while( candP[ left ].absRC > cut) + left++; + while( cut > candP[ right ].absRC) + right--; + + if( left < right ) + Swap( candP[ left ] , candP[ right ] ); + + if(left <= right) { + left++; + right--; + } + } while( left <= right ); + + if( min < right ) + SortPrimalCandidateList( min , right ); + if( ( left < max ) && ( left <= hotListSize ) ) + SortPrimalCandidateList( left , max ); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcDType* MCFSimplex::RuleDualCandidateListPivot( void ) +{ + Index next = 0; + // Check if the left arcs in the list continue to violate the primal condition + for( Index i = 2 ; ( i <= hotListSize ) && ( i <= tempCandidateListSize ) ; + i++ ) { + nodeDType *node = candD[ i ].node; + arcDType *arc = node->enteringTArc; + cFNumber flow = arc->flow; + if( LTZ( flow , EpsFlw ) ) { + next++; + candD[ next ].node = node; + candD[ next ].absInfeas = ABS( flow ); + } + + if( GT( flow , arc->upper , EpsFlw ) ) { + next++; + candD[ next ].node = node; + candD[ next ].absInfeas = flow - arc->upper; + } + } + + tempCandidateListSize = next; + Index oldGroupPos = groupPos; + // Search other arcs to fill the list + do { + nodeDType *node = nodesD + groupPos; + for( node ; node < stopNodesD ; node += numGroup ) { + arcDType *arc = node->enteringTArc; + cFNumber flow = arc->flow; + if( LTZ( flow , EpsFlw ) ) { + tempCandidateListSize++; + candD[ tempCandidateListSize ].node = node; + candD[ tempCandidateListSize ].absInfeas = ABS( flow ); + } + + if( GT( flow , arc->upper , EpsFlw) ) { + tempCandidateListSize++; + candD[ tempCandidateListSize ].node = node; + candD[ tempCandidateListSize ].absInfeas = flow - arc->upper; + } + } + + groupPos++; + if( groupPos == numGroup ) + groupPos = 0; + + } while( ( tempCandidateListSize < hotListSize ) && + ( groupPos != oldGroupPos ) ); + + if( tempCandidateListSize ) { + SortDualCandidateList( 1 , tempCandidateListSize ); + return( (candD[ 1 ].node)->enteringTArc ); + } + else + return( NULL ); + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFSimplex::InitializeDualCandidateList( void ) +{ + numGroup = ( ( n - 1 ) / numCandidateList ) + 1; + groupPos = 0; + tempCandidateListSize = 0; + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFSimplex::SortDualCandidateList(Index min, Index max) +{ + Index left = min; + Index right = max; + FNumber cut = candD[ ( left + right ) / 2 ].absInfeas; + do { + while( candD[ left ].absInfeas > cut ) + left++; + while( cut > candD[ right ].absInfeas ) + right--; + if( left < right ) + Swap( candD[left ] , candD[ right ] ); + if( left <= right) { + left++; + right--; + } + } while( left <= right ); + + if( min < right ) + SortDualCandidateList( min , right ); + if( (left < max) && ( left <= hotListSize ) ) + SortDualCandidateList( left , max ); + } + +/*--------------------------------------------------------------------------*/ + +template<class N, class RCT> +inline void MCFSimplex::AddPotential( N *r , RCT delta ) +{ + int level = r->subTreeLevel; + N *n = r; + + do { + n->potential = n->potential + delta; + n = n->nextInT; + } while ( ( n ) && ( n->subTreeLevel > level ) ); + } + +/*--------------------------------------------------------------------------*/ + +template<class N> +inline void MCFSimplex::ComputePotential( N *r ) +{ + N *n = r; + int level = r->subTreeLevel; + FONumber cost; + // If "n" is not the dummy root, the potential of "r" is computed. + // If "n" is the dummy root, the potential of dummy root is a constant. + + do { + if( n->enteringTArc ) { + cost = ( n->enteringTArc )->cost; + #if (QUADRATICCOST) + // Also field "sumQuadratic" is updated + n->sumQuadratic = ( Father( n , n->enteringTArc ) )->sumQuadratic + + ( n->enteringTArc )->quadraticCost; + + if( ! ETZ( ( n->enteringTArc )->flow , EpsFlw ) ) + cost = cost + ( ( n->enteringTArc )->quadraticCost * ( n->enteringTArc )->flow ); + #endif + + if( n == ( n->enteringTArc )->head ) + n->potential = ( Father( n , n->enteringTArc ) )->potential + cost; + else + n->potential = ( Father( n , n->enteringTArc ) )->potential - cost; + } + n = n->nextInT; + } while( ( n ) && ( n->subTreeLevel > level ) ); + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::CreateInitialPModifiedBalanceVector( void ) +{ + int i = 0; + delete[] modifiedBalance; + modifiedBalance = new FNumber[ n ]; + // Initialited every node's modifiedBalance to his balance + for ( nodePType *node = nodesP ; node != stopNodesP ; node++ ) { + modifiedBalance[i] = node->balance; + i++; + } + + // Modify the vector according to the arcs out of base with flow non zero + // Scan the real arcs + + for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) { + #if( QUADRATICCOST ) + if( ( ! ETZ( arc->flow , EpsFlw ) ) && + ( ( arc->tail )->enteringTArc != arc ) && + ( ( arc->head )->enteringTArc != arc ) ) { + i = (arc->tail) - nodesP; + modifiedBalance[ i ] += arc->flow; + i = (arc->head) - nodesP; + modifiedBalance[ i ] -= arc->flow; + } + #else + if( arc->ident == AT_UPPER ) { + i = (arc->tail) - nodesP; + modifiedBalance[ i ] += arc->upper; + i = (arc->head) - nodesP; + modifiedBalance[ i ] -= arc->upper; + } + #endif + } + + // Scan the dummy arcs + for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) { + #if( QUADRATICCOST ) + if ( ( ! ETZ( arc->flow , EpsFlw ) ) && + ( ( arc->tail )->enteringTArc != arc ) && + ( ( arc->head )->enteringTArc != arc ) ) { + i = (arc->tail) - nodesP; + modifiedBalance[ i ] += arc->flow; + i = (arc->head) - nodesP; + modifiedBalance[ i ] -= arc->flow; + } + #else + if (arc->ident == AT_UPPER) { + i = (arc->tail) - nodesP; + modifiedBalance[ i ] += arc->upper; + i = (arc->head) - nodesP; + modifiedBalance[ i ] -= arc->upper; + } + #endif + } + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::PostPVisit( nodePType *r ) +{ + // The method controls if "r" is a leaf in T + bool rLeaf = false; + int i = r - nodesP; + if( r->nextInT ) + if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel ) + rLeaf = true; + else + rLeaf = true; + + if( rLeaf ) // If "r" is a leaf + if( ( r->enteringTArc)->head == r ) // If enteringTArc of "r" goes in "r" + ( r->enteringTArc )->flow = modifiedBalance[ i ]; + else // If enteringTArc of "r" goes out "r" + ( r->enteringTArc )->flow = - modifiedBalance[ i ]; + else { // If "r" isn't a leaf + nodePType *desc = r->nextInT; + // Call PostPVisit for every child of "r" + while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) { + if( desc->subTreeLevel - 1 == r->subTreeLevel ) { // desc is a son of r + PostPVisit( desc ); + + if( ( desc->enteringTArc )->head == r ) // enteringTArc of desc goes in r + modifiedBalance[ i ] -= ( desc->enteringTArc )->flow; + else // If enteringTArc of "desc" goes out "r" + modifiedBalance[ i ] += ( desc->enteringTArc )->flow; + } + desc = desc->nextInT; + } + + if( r != dummyRootP ) + if( ( r->enteringTArc )->head == r ) // If enteringTArc of "r" goes in "r" + ( r->enteringTArc )->flow = modifiedBalance[ i ]; + else // If enteringTArc of "r" goes out "r" + ( r->enteringTArc )->flow = - modifiedBalance[ i ]; + } + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::BalanceFlow( nodePType *r ) +{ + // used only by Primal Simplex to restore a primal feasible solution. + if( r == dummyRootP ) { + nodePType *node = dummyRootP->nextInT; + while( node ) { + // call this function recursively for every son of dummy root + if( node->subTreeLevel == 1 ) + BalanceFlow( node ); + + node = node->nextInT; + } + } + else { + // The method controls if "r" is a leaf in T + bool rLeaf = false; + if( r->nextInT ) + if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel ) + rLeaf = true; + else + rLeaf = true; + + if( rLeaf ) // If "r" is a leaf + AdjustFlow( r ); // The method controls if entering basic arc in "r" is + // not feasible; in case adjust its flow + else { // If "r" isn't a leaf + nodePType *node = r->nextInT; + // Balance the flow of every child of "r" + while ( ( node ) && ( node->subTreeLevel > r->subTreeLevel ) ) { + if( node->subTreeLevel == r->subTreeLevel + 1 ) + BalanceFlow( node ); + node = node->nextInT; + } + + // The method controls if entering basic arc in "r" is not feasible; + //in case adjust its flow + AdjustFlow( r ); + } + } + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::AdjustFlow( nodePType *r ) +{ + arcPType *arc = r->enteringTArc; + if( arc >= dummyArcsP ) // If entering arc of "r" is a dummy arc + if( LTZ( arc->flow , EpsFlw ) ) { + // If this dummy arc has flow < 0, the algorithm overturns the arc + nodePType *temp = arc->tail; + arc->tail = arc->head; + arc->head = temp; + arc->flow = -arc->flow; + } + else { // If entering arc of "r" is not a dummy arc + bool orientationDown = ( arc->head == r ); + FNumber delta = 0; + if( LTZ( arc->flow , EpsFlw ) ) { // If flow is < 0 + delta = -arc->flow; + arc->flow = 0; + #if( ! QUADRATICCOST ) + arc->ident = AT_LOWER; + #endif + } + + if( GT( arc->flow , arc->upper , EpsFlw ) ) { + // If flow goes over the capacity of the arc + delta = arc->upper - arc->flow; + arc->flow = arc->upper; + #if( ! QUADRATICCOST ) + arc->ident = AT_UPPER; + #endif + } + + /* This arc goes out from the basis, and the relative dummy arc goes in T. + Then the algorithm push flow in the cycle made by the arc and some arcs + of T to balance the flow. */ + + if( ! ETZ( delta , EpsFlw ) ) { + nodePType *node = Father( r , arc ); + while( node != dummyRootP ) { + arc = node->enteringTArc; + if( ( arc->head == node ) == orientationDown ) + arc->flow += delta; + else + arc->flow -= delta; + + node = Father( node , arc ); + } + + arcPType *dummy = dummyArcsP + ( r - nodesP ); + #if( ! QUADRATICCOST ) + dummy->ident = BASIC; + #endif + + /* Update the structure of the tree. If entering basic arc of "r" is + changed, subtree of "r"is moved next dummy root. */ + r->enteringTArc = dummy; + int deltaLevel = 1 - r->subTreeLevel; + nodePType *lastNode = CutAndUpdateSubtree( r , deltaLevel ); + PasteSubtree( r , lastNode , dummyRootP ); + if( ( dummy->head == r ) != orientationDown ) + dummy->flow += delta; + else + dummy->flow -= delta; + + if( LTZ( dummy->flow , EpsFlw ) ) { + nodePType *temp = dummy->tail; + dummy->tail = dummy->head; + dummy->head = temp; + dummy->flow = -dummy->flow; + } + } + } + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::CreateInitialDModifiedBalanceVector( void ) +{ + #if( ! QUADRATICCOST ) + int i = 0; + modifiedBalance = new FNumber[ n ]; + // Initialited every node's modifiedBalance to his balance + for( nodeDType *node = nodesD ; node != stopNodesD ; node++ ) { + modifiedBalance[ i ] = node->balance; + i++; + } + + // Modify the vector according to the arcs out of base with flow non zero + // Scan the real arcs + for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) + if( arc->ident == AT_UPPER ) { + i = (arc->tail) - nodesD; + modifiedBalance[ i ] += arc->upper; + i = (arc->head) - nodesD; + modifiedBalance[ i ] -= arc->upper; + } + + // Scan the dummy arcs + for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ ) + if( arc->ident == AT_UPPER ) { + i = (arc->tail) - nodesD; + modifiedBalance[ i ] += arc->upper; + i = (arc->head) - nodesD; + modifiedBalance[ i ] -= arc->upper; + } + #endif + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::PostDVisit( nodeDType *r ) +{ + #if( ! QUADRATICCOST ) + // The method controls if "r" is a leaf in T + bool rLeaf = false; + int i = r - nodesD; + if( r->nextInT ) + if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel ) + rLeaf = true; + else + rLeaf = true; + + if( rLeaf ) // If "r" is a leaf + if( ( r->enteringTArc)->head == r ) // If enteringTArc of "r" goes in "r" + ( r->enteringTArc )->flow = modifiedBalance[ i ]; + else // If enteringTArc of "r" goes out "r" + ( r->enteringTArc )->flow = - modifiedBalance[ i ]; + else { // If "r" isn't a leaf + nodeDType *desc = r->nextInT; + // Call PostDVisit for every child of "r" + while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) { + if( desc->subTreeLevel -1 == r->subTreeLevel ) { // desc is a son of r + PostDVisit( desc ); + + if( ( desc->enteringTArc )->head == r ) // enteringTArc of desc goes in r + modifiedBalance[ i ] -= ( desc->enteringTArc )->flow; + else // If enteringTArc of "desc" goes out "r" + modifiedBalance[ i ] += ( desc->enteringTArc )->flow; + } + + desc = desc->nextInT; + } + + if( r != dummyRootD ) + if( ( r->enteringTArc )->head == r ) // If enteringTArc of "r" goes in "r" + ( r->enteringTArc )->flow = modifiedBalance[ i ]; + else // If enteringTArc of "r" goes out "r" + ( r->enteringTArc )->flow = - modifiedBalance[ i ]; + } + #endif + } + +/*--------------------------------------------------------------------------*/ + +inline void MCFSimplex::ResetWhenInT2( void ) +{ + for( nodeDType *n = nodesD ; n != stopNodesD ; n++) + n->whenInT2 = 0; + } + +/*--------------------------------------------------------------------------*/ + +template<class N, class A> +inline N* MCFSimplex::Father( N *n , A *a ) +{ + if( a == NULL ) + return NULL; + + if( a->tail == n ) + return( a->head ); + else + return( a->tail ); + } + +/*-------------------------------------------------------------------------*/ + +inline MCFSimplex::FONumber MCFSimplex::GetFO( void ) +{ + FONumber fo = 0; + if( usePrimalSimplex ) { + arcPType *arco; + for( arco = arcsP ; arco != stopArcsP ; arco++ ) { + #if( QUADRATICCOST ) + if( ! ETZ( arco->flow , EpsFlw ) ) + fo += arco->flow * ( arco->cost + arco->flow * arco->quadraticCost / 2 ); + #else + if( ( arco->ident == BASIC ) || ( arco->ident == AT_UPPER ) ) + fo += arco->cost * arco->flow; + #endif + } + + for( arco = dummyArcsP ; arco != stopDummyP ; arco++ ) { + #if( QUADRATICCOST ) + if( ! ETZ( arco->flow , EpsFlw ) ) + fo += arco->flow * ( arco->cost + arco->flow * arco->quadraticCost / 2 ); + #else + if( ( arco->ident == BASIC ) || ( arco->ident == AT_UPPER ) ) + fo += arco->cost * arco->flow; + #endif + } + } + else { + arcDType *a; + for( a = arcsD ; a != stopArcsD ; a++ ) { + #if (QUADRATICCOST) + fo += ( a->cost * a->flow ) + ( a->quadraticCost * a->flow * a->flow ) / 2; + #else + if( ( a->ident == BASIC ) || (a->ident == AT_UPPER ) ) + fo += a->cost * a->flow; + #endif + } + + for( a = dummyArcsD ; a != stopDummyD ; a++) { + #if (QUADRATICCOST) + fo += ( a->cost * a->flow ) + ( a->quadraticCost * a->flow * a->flow ) / 2; + #else + if( ( a->ident == BASIC ) || ( a->ident == AT_UPPER ) ) + fo += a->cost * a->flow; + #endif + } + } + + return( fo ); + } + +/*-------------------------------------------------------------------------*/ + +void MCFSimplex::PrintPNode( nodePType *nodo ) +{ + if( nodo ) + if( nodo != dummyRootP ) + cout << ( nodo - nodesP + 1 ); + else + cout << "r"; + else + cout << ".."; + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::PrintPArc( arcPType *arc ) +{ + if( arc ) { + cout << "("; + PrintPNode( arc->tail ); + cout << ", "; + PrintPNode( arc->head ); + cout << ")"; + } + else + cout << ".."; +} + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::PrintDNode( nodeDType *nodo ) +{ + if( nodo ) + if( nodo != dummyRootD ) + cout << ( nodo - nodesD + 1 ); + else + cout << "r"; + else + cout << ".."; + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::PrintDArc( arcDType *arc ) +{ + if( arc ) { + cout << "("; + PrintDNode( arc->tail ); + cout << ", "; + PrintDNode( arc->head ); + cout << ")"; + } + else + cout << ".."; + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::nodePType* MCFSimplex::RecoverPNode( Index ind ) +{ + if( ( ind < 0 ) || ( ind > n ) ) + return( NULL ); + if( ind ) + return( nodesP + ind - 1 ); + else + return( dummyRootP ); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcPType* MCFSimplex::RecoverPArc( nodePType *tail , + nodePType *head ) +{ + if( ( tail == NULL ) || ( head == NULL ) ) + return( NULL ); + + arcPType *arc = arcsP; + while( ( arc->tail != tail ) || ( arc->head != head ) ) { + arc++; + if( arc == stopArcsP ) + arc = dummyArcsP; + if( arc == stopDummyP ) + return( NULL ); + } + + return( arc ); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::nodeDType* MCFSimplex::RecoverDNode( Index ind ) +{ + if( ( ind < 0 ) || ( ind > n ) ) + return( NULL ); + + if( ind ) + return( nodesD + ind - 1 ); + else + return( dummyRootD ); + } + +/*--------------------------------------------------------------------------*/ + +MCFSimplex::arcDType* MCFSimplex::RecoverDArc( nodeDType *tail , + nodeDType *head ) +{ + if( ( tail == NULL ) || ( head == NULL ) ) + return( NULL ); + + arcDType *arc = arcsD; + while( ( arc->tail != tail ) || ( arc->head != head ) ) { + arc++; + if( arc == stopArcsD ) + arc = dummyArcsD; + if( arc == stopDummyD ) + return( NULL ); + } + + return( arc ); + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::infoPNode( nodePType *node , int tab ) +{ + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + cout << "Nodo "; + PrintPNode( node ); + cout << ": b = " << node->balance << " y = " << node->potential << endl; + #if( UNIPI_VIS_NODE_BASIC_ARC ) + cout << ": TArc="; + PrintPArc( node->enteringTArc ); + cout << endl; + #endif + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::infoPArc( arcPType *arc , int ind , int tab ) +{ + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + cout << "Arco "; + PrintPArc( arc ); + cout << ": x = " << arc->flow; + #if( UNIPI_VIS_ARC_UPPER ) + cout << " u = " << arc->upper; + #endif + #if( UNIPI_VIS_ARC_COST ) + cout << " c = " << arc->cost; + #endif + #if( QUADRATICCOST ) + #if( UNIPI_VIS_ARC_Q_COST ) + cout << " q = " << arc->quadraticCost; + #endif + cout << endl; + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + #if( UNIPI_VIS_ARC_REDUCT_COST ) + cout << " rc = " << MCFGetRC( ind ); + #endif + #else + cout << endl; + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + #if( UNIPI_VIS_ARC_REDUCT_COST ) + cout << " rc = " << MCFGetRC( ind ); + #endif + #if( UNIPI_VIS_ARC_STATE ) + switch( arc->ident ) { + case( BASIC ): cout << " in T"; break; + case( AT_LOWER ): cout << " in L"; break; + case( AT_UPPER ): cout << " in U"; break; + case( DELETED ): cout << " canceled"; break; + case( CLOSED ): cout << " closed"; + } + #endif + #endif + cout << endl; + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::infoDNode( nodeDType *node , int tab ) +{ + for( int t = 0 ; t < tab; t++ ) + cout << "\t"; + cout << "Nodo "; + PrintDNode( node ); + cout << ": b = " << node->balance << " y = " << node->potential; + #if( UNIPI_VIS_NODE_BASIC_ARC ) + cout << ": TArc="; + PrintDArc( node->enteringTArc ); + cout << endl; + #endif + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::infoDArc( arcDType *arc , int ind , int tab ) +{ + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + cout << "Arco "; + PrintDArc( arc ); + cout << " x = " << arc->flow; + #if( UNIPI_VIS_ARC_UPPER ) + cout << " u = " << arc->upper; + #endif + #if( UNIPI_VIS_ARC_COST ) + cout << " c = " << arc->cost; + #endif + #if( QUADRATICCOST ) + #if( UNIPI_VIS_ARC_Q_COST ) + cout << " q = " << arc->quadraticCost; + #endif + cout << endl; + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + #if( UNIPI_VIS_ARC_REDUCT_COST ) + cout << " rc = " << MCFGetRC( ind ); + #endif + #else + cout << endl; + for( int t = 0 ; t < tab ; t++ ) + cout << "\t"; + #if( UNIPI_VIS_ARC_REDUCT_COST ) + cout << " rc = " << MCFGetRC( ind ); + #endif + #if (UNIPI_VIS_ARC_STATE) + switch( arc->ident ) { + case( BASIC ): cout << " in T"; break; + case( AT_LOWER ): cout << " in L"; break; + case( AT_UPPER ): cout << " in U"; break; + case( DELETED ): cout << " canceled"; break; + case( CLOSED ): cout << " closed"; + } + #endif + #endif + cout << endl; + } + +/*--------------------------------------------------------------------------*/ + +void MCFSimplex::ShowSituation( int tab ) +{ + if( usePrimalSimplex ) { + arcPType *arc; + nodePType *node; + int i = 0; + for( arc = arcsP ; arc != stopArcsP ; arc++ ) { + infoPArc( arc , i , tab ); + i++; + } + cout << endl; + #if( UNIPI_VIS_DUMMY_ARCS ) + i = 0; + for( arc = dummyArcsP ; arc != stopDummyP ; arc++ ) { + infoPArc( arc , i , tab ); + i++; + } + cout << endl; + #endif + infoPNode( dummyRootP , tab ); + for( node = nodesP ; node != stopNodesP ; node++ ) + infoPNode( node , tab ); + } + else { + arcDType *arc; + nodeDType *node; + int i = 0; + for( arc = arcsD ; arc != stopArcsD ; arc++ ) { + infoDArc( arc , i , tab ); + i++; + } + cout << endl; + #if( UNIPI_VIS_DUMMY_ARCS ) + i = 0; + for( arc = dummyArcsD ; arc != stopDummyD ; arc++) { + infoDArc( arc , i , tab ); + i++; + } + cout << endl; + #endif + infoDNode( dummyRootD , tab ); + for( node = nodesD ; node != stopNodesD ; node++ ) + infoDNode( node , tab ); + } + } + +/*-------------------------------------------------------------------------*/ +/*---------------------- End File MCFSimplex.C ----------------------------*/ +/*-------------------------------------------------------------------------*/ diff --git a/impl/MCFSimplex.h b/impl/MCFSimplex.h new file mode 100644 index 0000000..a8d84f6 --- /dev/null +++ b/impl/MCFSimplex.h @@ -0,0 +1,1086 @@ +/*--------------------------------------------------------------------------*/ +/*------------------------- File MCFSimplex.h ------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @file + * Linear and Quadratic Min Cost Flow problems solver based on the (primal and + * dual) simplex algorithm. Conforms to the standard MCF interface defined in + * MCFClass.h. + * + * \Version 1.00 + * + * \date 29 - 08 - 2011 + * + * \author Alessandro Bertolini \n + * Antonio Frangioni \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * Copyright © 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni + */ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*----------------------------- DEFINITIONS --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef __MCFSimplex + #define __MCFSimplex /* self-identification: #endif at the end of the file */ + +/*--------------------------------------------------------------------------*/ +/*------------------------------ INCLUDES ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#include "MCFClass.h" + +/*@} -----------------------------------------------------------------------*/ +/*--------------------------- NAMESPACE ------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) +namespace MCFClass_di_unipi_it +{ +#endif + +/*--------------------------------------------------------------------------*/ +/*-------------------------------- MACROS ----------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup MCFSimplex_MACROS Compile-time switches in MCFSimplex.h + There is only one macro in MCFSimplex, but it is very important! + @{ */ + +#define QUADRATICCOST 0 + +/**< Setting QUADRATICCOST == 1 the solver can solve problems with linear and + quadratic costs too (but the latter only with the Primal Simplex). + The reason for having a macro is that when quadratic costs are present the + "arcType" struct has the additional field "quadraticCost" to hold it. + Furthermore, the field "ident" is not created because the solver doesn't + use the classical TLU tripartition. Instead, closed arcs and deleted arcs + are characterized as follows: + - closed arcs have the field "cost" to INFINITY (Inf<FNumber>()); + - deleted arcs have the field "upper" to INFINITY and the "tail" and + "head" field are NULL. + Furthermore, the solver needs the variables "ignoredEnteringArc" and + "firstIgnoredEnteringArc", used to avoid nasty loops during the execution + of the Quadratic Primal Simplex algorithm. + If, instead, QUADRATICCOST == 0 then the solver can solve only problems + with linear costs. Hence, the field "quadraticCost" is useless and it + isn't created. Furthermore, Primal Simplex and Dual Simplex use the + tripartition TLU to divide the arcs, so the solver creates the field + "ident", which differentiates the set of the arcs in among deleted arcs, + closed arcs, arcs in T, arcs in L, arcs in U. + Thus, with QUADRATICCOST == 0 the solver cannot solve problems with + quadratic costs, but it does solve problems with linear costs faster. */ + +/*@} end( group( MCFCLASS_MACROS ) ) */ + +/*--------------------------------------------------------------------------*/ +/*---------------------------- CLASSES -------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup MCFSimplex_CLASSES Classes in MCFSimplex.h + @{ */ + +/** The MCFSimplex class derives from the abstract base class MCFClass, thus + sharing its (standard) interface, and implements both the Primal and + Dual network simplex algorithms for solving (Linear and Quadratic) + Min Cost Flow problems */ + +class MCFSimplex: public MCFClass +{ + +/*--------------------------------------------------------------------------*/ +/*----------------------- PUBLIC PART OF THE CLASS -------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-- --*/ +/*-- The following methods and data are the actual interface of the --*/ +/*-- class: the standard user should use these methods and data only. --*/ +/*-- --*/ +/*--------------------------------------------------------------------------*/ + + public: + +/*--------------------------------------------------------------------------*/ +/*---------------------------- PUBLIC TYPES --------------------------------*/ +/*--------------------------------------------------------------------------*/ + +/** Public enum describing the parameters of MCFSimplex. */ + + enum SimplexParam + { + kAlgPrimal = kLastParam , ///< parameter to set algorithm (Primal/Dual): + kAlgPricing , ///< parameter to set algorithm of pricing + kNumCandList , /**< parameter to set the number of candidate + list for Candidate List Pivot method */ + kHotListSize , /**< parameter to set the size of Hot List + for Candidate List Pivot method */ + kRecomputeFOLimits , /**< parameter to set the number of iterations + in which quadratic Primal Simplex computes + "manually" the f.o. value */ + kEpsOpt /**< parameter to set the precision of the + objective function value for the + quadratic Primal Simplex */ + }; + +/** Public enum describing the pricing rules in MCFSimplex::SetAlg(). */ + + enum enumPrcngRl + { + kDantzig = 0, ///< Dantzig's rule (most violated constraint) + kFirstEligibleArc , ///< First eligible arc in round-robin + kCandidateListPivot ///< Candidate List Pivot Rule + }; + +/*--------------------------------------------------------------------------*/ +/*--------------------------- PUBLIC METHODS -------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*---------------------------- CONSTRUCTOR ---------------------------------*/ +/*--------------------------------------------------------------------------*/ + + MCFSimplex( cIndex nmx = 0 , cIndex mmx = 0 ); + +/**< Constructor of the class, as in MCFClass::MCFClass(). */ + +/*--------------------------------------------------------------------------*/ +/*-------------------------- OTHER INITIALIZATIONS -------------------------*/ +/*--------------------------------------------------------------------------*/ + + void LoadNet( cIndex nmx = 0 , cIndex mmx = 0 , cIndex pn = 0 , + cIndex pm = 0 , cFRow pU = NULL , cCRow pC = NULL , + cFRow pDfct = NULL , cIndex_Set pSn = NULL , + cIndex_Set pEn = NULL ); + +/**< Inputs a new network, as in MCFClass::LoadNet(). */ + +/*--------------------------------------------------------------------------*/ + + void SetAlg( bool UsPrml , char WhchPrc ); + +/**< Selects which algorithm (Primal vs Dual Network Simplex), and which + pricing rule within the algorithm, is used. + + If UsPrml == TRUE then the Primal Network Simplex algorithm is used, + otherwise the Dual Network Simplex is used. + + The allowed values for WhchPrc are: + + - kDantzig Dantzig's pricing rule, i.e., most violated dual + constraint; this can only be used with the Primal + Network Simplex + + - kFirstEligibleArcA "dumb" rule, first eligible arc in round-robin; + + - kCandidateListPivot Candidate List Pivot Rule + + If this method is *not* called, the Primal Network Simplex with the + Candidate List Pivot Rule (the best setting on most instances) is + used. */ + +/*--------------------------------------------------------------------------*/ + + void SetPar( int par , int val ); + +/**< Set general integer parameters. + + @param par is the parameter to set; since this method accepts an int + value, the enum SimplexParam can be used in addition to the + enum MCFParam to specify the integer parameters (every enum + is an int). + + @param value is the value to assign to the parameter. */ + +/*-------------------------------------------------------------------------*/ + + void SetPar( int par , double val ); + +/**< Set general float parameters. + + @param par is the parameter to set; since this method accepts an int + value, the enum SimplexParam can be used in addition to the + enum MCFParam to specify the float parameters (every enum + is an int). + + @param value is the value to assign to the parameter. */ + +/*--------------------------------------------------------------------------*/ +/*-------------------- METHODS FOR SOLVING THE PROBLEM ---------------------*/ +/*--------------------------------------------------------------------------*/ + + void SolveMCF( void ); + +/*--------------------------------------------------------------------------*/ +/*---------------------- METHODS FOR READING RESULTS -----------------------*/ +/*--------------------------------------------------------------------------*/ + + void MCFGetX( FRow F , Index_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + +/*--------------------------------------------------------------------------*/ + + void MCFGetRC( CRow CR , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ) ; + + CNumber MCFGetRC( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFGetPi( CRow P , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + +/*--------------------------------------------------------------------------*/ + + FONumber MCFGetFO( void ); + +/*--------------------------------------------------------------------------*/ +/*-------------- METHODS FOR READING THE DATA OF THE PROBLEM ---------------*/ +/*--------------------------------------------------------------------------*/ + + virtual void MCFArcs( Index_Set Startv , Index_Set Endv , + cIndex_Set nms = NULL , cIndex strt = 0 , + Index stp = Inf<Index>() ); + + inline Index MCFSNde( cIndex i ); + + inline Index MCFENde( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFCosts( CRow Costv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + inline CNumber MCFCost( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFQCoef( CRow Qv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + inline CNumber MCFQCoef( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + inline FNumber MCFUCap( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + inline FNumber MCFDfct( cIndex i ); + +/*--------------------------------------------------------------------------*/ +/*------------- METHODS FOR ADDING / REMOVING / CHANGING DATA --------------*/ +/*--------------------------------------------------------------------------*/ +/*------- Changing the costs, QCoef, deficits and upper capacities ---------*/ +/*--------------------------------------------------------------------------*/ + + void ChgCosts( cCRow NCost , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + void ChgCost( Index arc , cCNumber NCost ); + +/*--------------------------------------------------------------------------*/ + + void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + void ChgQCoef( Index arc , cCNumber NQCoef ); + +/*--------------------------------------------------------------------------*/ + + void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + void ChgDfct( Index nod , cFNumber NDfct ); + +/*--------------------------------------------------------------------------*/ + + void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf<Index>() ); + + void ChgUCap( Index arc , cFNumber NCap ); + +/*--------------------------------------------------------------------------*/ +/*--------------- Modifying the structure of the graph ---------------------*/ +/*--------------------------------------------------------------------------*/ + + void CloseArc( cIndex name ); + + void DelNode( cIndex name ); + + bool IsClosedArc( cIndex name ); + + void OpenArc( cIndex name ) ; + + Index AddNode( cFNumber aDfct ); + + void ChangeArc( cIndex name , cIndex nSS = Inf<Index>() , + cIndex nEN = Inf<Index>() ); + + void DelArc( cIndex name ); + + Index AddArc( cIndex Start , cIndex End , cFNumber aU , cCNumber aC ); + + bool IsDeletedArc( cIndex name ); + +/*--------------------------------------------------------------------------*/ +/*------------------------------ DESTRUCTOR --------------------------------*/ +/*--------------------------------------------------------------------------*/ + + ~MCFSimplex(); + +/*--------------------------------------------------------------------------*/ +/*--------------------- PRIVATE PART OF THE CLASS --------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-- --*/ +/*-- Nobody should ever look at this part: everything that is under this --*/ +/*-- advice may be changed without notice in any new release of the code. --*/ +/*-- --*/ +/*--------------------------------------------------------------------------*/ + + private: + +/*--------------------------------------------------------------------------*/ +/*-------------------------- PRIVATE DATA TYPES ----------------------------*/ +/*--------------------------------------------------------------------------*/ +/*-- --*/ +/*-- Let T \subseteq A be a spanning tree, and consider some node v \in V --*/ +/*-- \setminus { 0 }. There is an unique (undirected) path, denoted by --*/ +/*-- P(v), defined by T from v to the root node 0. The arc in P(v), which --*/ +/*-- is incident to v, is called the *basic arc* of v. The other terminal --*/ +/*-- node u of this basic arc is called the *father node* of v. The --*/ +/*-- spanning tree T is represented saving the basic arc of every node, --*/ +/*-- and maintaining the order of the nodes and the depth as to T root --*/ +/*-- after a Post-Visit of T. This order is saved in a bidirectional list --*/ +/*-- written in the node. --*/ +/*-- --*/ +/*-- The Primal Simplex uses a different data structure than the Dual --*/ +/*-- Simplex, because the Dual Simplex needs additional data (mainly the --*/ +/*-- Backward Star and Forward Star. Furthermore, the Primal Simplex uses --*/ +/*-- different data structures in the quadratic case. --*/ +/*-- --*/ +/*--------------------------------------------------------------------------*/ + + struct arcPType; // pre-declaration of the arc structure (pointers to arcs + // are contained in the node structure) for Primal Simplex + + struct arcDType; // pre-declaration of the arc structure (pointers to arcs + // are contained in the node structure) for Dual Simplex + + typedef double iteratorType; // type for the iteration counter and array + // "whenInT2" + + struct nodePType { // node structure for Primal Simplex - - - - - - - - + nodePType *prevInT; // previous node in the order of the Post-Visit on T + + nodePType *nextInT; // next node in the order of the Post-Visit on T + + arcPType *enteringTArc; // entering basic arc of this node + + FNumber balance; // supply/demand of this node; a node is called a + // supply node, a demand node, or a transshipment + // node depending upon whether balance is larger + // than, smaller than, or equal to zero + #if( QUADRATICCOST ) + CNumber sumQuadratic; // the sum of the quadratic coefficients of the tree's arcs + // from root of T to the node + + FONumber potential; // the node potential corresponding with the flow + // conservation constrait of this node + #else + CNumber potential; // the node potential corresponding with the flow + // conservation constrait of this node + #endif + + int subTreeLevel; // the depth of the node in T as to T root + + }; // end( struct( nodePType ) ) + + struct nodeDType { // node structure for Dual Simplex - - - - - - - - - + nodeDType *prevInT; // previous node in the order of the Post-Visit on T + + nodeDType *nextInT; // next node in the order of the Post-Visit on T + + arcDType *enteringTArc; // entering basic arc of this node + + FNumber balance; // supply/demand of this node; a node is called a + // supply node, a demand node, or a transshipment + // node depending upon whether balance is larger + // than, smaller than, or equal to zero + + #if( QUADRATICCOST ) + CNumber sumQuadratic; // the sum of the quadratic coefficients of the tree's arcs + // from root of T to the node + + FONumber potential; // the node potential corresponding with the flow + // conservation constrait of this node + #else + CNumber potential; // the node potential corresponding with the flow + // conservation constrait of this node + #endif + + int subTreeLevel; // the depth of the node in T as to T root + iteratorType whenInT2; // the last iteration where a node is in subtree T2 + + Index numArcs; // the number of the arcs which enter/exit from node + arcDType *firstBs; // the first arc in the node's Backward Star + arcDType *firstFs; // the first arc in the node's Forward Star + + }; // end( struct( nodeDType ) ) + + struct arcPType { // arc structure for Primal Simplex - - - - - - - - + nodePType *tail; // tail node + nodePType *head; // head node + FNumber flow; // arc flow + CNumber cost; // arc linear cost + + #if( QUADRATICCOST ) + CNumber quadraticCost; // arc quadratic cost + #else + char ident; // tells if arc is deleted, closed, in T, L, or U + #endif + + FNumber upper; // arc upper bound + }; // end( struct( arcPType ) ) + + struct arcDType { // arc structure for Primal Simplex - - - - - - - - + nodeDType *tail; // tail node + nodeDType *head; // head node + FNumber flow; // arc flow + CNumber cost; // arc linear cost + + #if( QUADRATICCOST ) + CNumber quadraticCost; // arc quadratic cost + #else + char ident; // indicates if arc is deleted, closed, in T, in L, or in U + #endif + + FNumber upper; // arc upper bound + arcDType *nextBs; // the next arc in the Backward Star of the arc's head + arcDType *nextFs; // the next arc in the Forward Star of the arc's tail + + }; // end( struct( arcDType ) ) + + struct primalCandidType { // Primal Candidate List- - - - - - - - - - - - - + arcPType *arc; // pointer to the violating primal bound arc + + #if(QUADRATICCOST ) + FONumber absRC; // absolute value of the arc's reduced cost + #else + CNumber absRC; // absolute value of the arc's reduced cost + #endif + }; // end( struct( primalCandidateType ) ) + + struct dualCandidType { // Dual Candidate List- - - - - - - - - - - - - - + nodeDType *node; // deepest node violating the dual bound arc + FNumber absInfeas; // absolute value of the arc's flow infeasibility + }; // end( struct( dualCandidateType ) ) + +/*--------------------------------------------------------------------------*/ +/*----------------------- PRIVATE DATA STRUCTURES -------------------------*/ +/*--------------------------------------------------------------------------*/ + + bool usePrimalSimplex; // TRUE if the Primal Network Simplex is used + + char pricingRule; // which pricing rule is used + + nodePType *nodesP; // vector of nodes: points to the n + 1 node + // structs (including the dummy root node) + // where the first node is indexed by zero + // and the last node is the dummy root node + + nodePType *dummyRootP; // the dummy root node + + nodePType *stopNodesP; // first infeasible node address = nodes + n + + arcPType *arcsP; // vector of arcs; this variable points to + // the m arc structs. + + arcPType *dummyArcsP; // vector of artificial dummy arcs: points + // to the artificial dummy arc variables and + // contains n arc structs. The artificial + // arcs are used to build artificial feasible + // starting bases. For each node i, there is + // exactly one dummy arc defined to connect + // the node i with the dummy root node. + + arcPType *stopArcsP; // first infeasible arc address = arcs + m + + arcPType *stopDummyP; // first infeasible dummy arc address + // = arcs + m + n + + arcPType *arcToStartP; // Dantzig Rule and First Eligible Arc Rule + // start their search from this arc + + nodeDType *nodesD; // vector of nodes: points to the n + 1 node + // structs (including the dummy root node) + // where the first node is indexed by zero + // and the last node is the dummy root node + + nodeDType *dummyRootD; // the dummy root node + + nodeDType *stopNodesD; // first infeasible node address = nodes + n + + arcDType *arcsD; // vector of arcs; this variable points to + // the m arc structs. + + arcDType *dummyArcsD; // vector of artificial dummy arcs: points to + // to the artificial dummy arc variables and + // contains n arc structs. The artificial + // arcs are used to build artificial feasible + // starting bases. For each node i, there is + // exactly one dummy arc defined to connect + // the node i with the dummy root node. + + arcDType *stopArcsD; // first infeasible arc address = arcs + m + + arcDType *stopDummyD; // first infeasible dummy arc address + // = arcs + m + n + + arcDType *arcToStartD; // Dantzig Rule and First Eligible Arc Rule + // start their search from this arc + + iteratorType iterator; // the current number of iterations + + primalCandidType *candP; // every element points to an element of the + // arcs vector which contains an arc violating + // dual bound + + dualCandidType *candD; // every element points to an element of the + // arcs vector which contains an arc violating + // primal bond + + Index numGroup; // number of the candidate lists + + Index tempCandidateListSize; // hot list dimension (it is variable) + + Index groupPos; // contains the actual candidate list + + Index numCandidateList; // number of candidate lists + + Index hotListSize; // number of candidate lists and hot list dimension + + Index forcedNumCandidateList; // value used to force the number of candidate list + + Index forcedHotListSize; // value used to force the number of candidate list + // and hot list dimension + + bool newSession; // true if algorithm is just started + + CNumber MAX_ART_COST; // large cost for artificial arcs + + FNumber *modifiedBalance; // vector of balance used by the PostVisit + + FONumber EpsOpt; // the precision of the objective function value + // for the quadratic case of the Primal Simplex + + int recomputeFOLimits; // after how many iterations the quadratic Primal + // Simplex computes "manually" the o.f. value + + #if( QUADRATICCOST ) + FONumber foValue; // the temporary objective function value + #endif + +/*--------------------------------------------------------------------------*/ +/*-------------------------- PRIVATE METHODS -------------------------------*/ +/*--------------------------------------------------------------------------*/ + + void MemAlloc( void ); + +/**< Method to allocate memory for the main data structures. It creates the + vector of nodes, the vector of arcs (real and dummy). If the algorithm is + using the Dual Simplex, it creates also the vectors whenInT2, firstIn and + nextIn, usefull to identify the next entering. */ + +/*--------------------------------------------------------------------------*/ + + void MemDeAlloc( bool whatDeAlloc ); + +/**< Method to deallocate memory for the main data structures created in + MemAlloc(). */ + +/*--------------------------------------------------------------------------*/ + + void MemAllocCandidateList( void ); + +/**< Method to allocate memory for the data structure used by the Candidate + List Pivot Rule. It creates the vector of candP (or candD in the Dual + Simplex case), determining the size of this vector on the basis of number + of arcs or nodes. */ + +/*--------------------------------------------------------------------------*/ + + void MemDeAllocCandidateList( void ); + +/**< Method to deallocate memory for the data structures created in + MemAllocCandidateList(). */ + +/*--------------------------------------------------------------------------*/ + + void CreateInitialPrimalBase( void ); + +/**< Method to create an initial feasible primal base. Add one node (dummyRoot) + to the network and connect this dummy root to each real node with n dummy + arcs. For each source node i, a dummy arc (i, r) exists. For each sink or + transit node j, a dummy arc (r, j) exists. The source nodes send their flows + to the dummy root through their dummy arcs, so the dummy root send all to + the sink nodes. The dummy root balance is 0, the costs of dummy arcs are + fixed to "infinity". */ + +/*--------------------------------------------------------------------------*/ + + void CreateInitialDualBase( void ); + +/**< Method to create an initial feasible dual base. Add one node (dummyRoot) + to the network and connect this dummy root to each real node with n dummy + arcs. The source nodes send their flows to the dummy root through their + dummy arcs, so the dummy root send all to the sink nodes. The dummy root + balance is 0, the costs of dummy arcs are fixed to "infinity". */ + +/*--------------------------------------------------------------------------*/ + + void CreateAdditionalDualStructures( void ); + +/**< The Dual Simplex needs nodes' Backward and Forward Star to work. So when + the Primal Simplex runs, these structures don't exist. When the Dual + Simplex starts, these structure are created in this method. */ + +/*--------------------------------------------------------------------------*/ + + void PrimalSimplex( void ); + +/**< Main method to implement the Primal Simplex algorithm. */ + +/*--------------------------------------------------------------------------*/ + + void DualSimplex( void ); + +/**< Main method to implement the Dual Simplex algorithm. */ + +/*--------------------------------------------------------------------------*/ + + template<class N, class A> + void UpdateT( A *h , A *k , N *h1 , N *h2 , N *k1 , N *k2 ); + + /**< Method to update the spanning tree T every iteration. + The spanning tree T is implemented with a bidirectional list stored + in the node structure, which represents the nodes' order after a + Post-Visit. The parameter of the method are the outgoing arc "h", the + incoming arc "k", and the four nodes on the extremity of these arcs + (for example h2 is the deepest node of the outgoing arc). Removing the + arc "h" splits T in two subtrees: T1 (which contains the root of T) and + T2, which will be re-connected by the incoming arc "k". + T2 will be reordered; in fact, the node "k2" becomes the root of T2 + instead of "h2" and the hierarchy of T2 will be overturned. Then T2 will + be moved; the root of T2 is changed, therefore the predecessor of the + root will become the node "k1". + + This method uses the methods cutAndUpdateSubtree() and pasteSubtree(). + First it cuts the node "k2" and its subtree from T using + cutAndUpdateSubtree(). Moreover the method cutAndUpdateSubtree() + updates the field "subTreeLevel" of every subtree's nodes, since k2's + subtree will be moved from the bottom to the top of T2. Then the method + pasteSubtree() puts this subtree in the bidirectional list after the node + "k1". The same operations will be applied to the old precedessor of "k2" + (which will become one of the childs of "k2"). This second subtree will + be cut, the subTreeLevel fields will be updated, and it will be inserted + in the bidirectional list after the k2's subtree. This is iterated until + the node "h2" is reached. */ + +/*--------------------------------------------------------------------------*/ + + template<class N> + N* CutAndUpdateSubtree( N *root, int delta ); + +/**< This method cuts a generic subtree from the spanning tree T. Then it + updates the field "subTreeLevel" of every subtree's nodes adding the value + "delta". This method returns the last node of the subtree. */ + +/*--------------------------------------------------------------------------*/ + + template<class N> + void PasteSubtree( N *root , N *lastNode , N *previousNode ); + +/**< This method inserts a generic subtree with root passed by parameter into + the spanning tree T, between the nodes "previousNode" and "lastNode". */ + +/*--------------------------------------------------------------------------*/ + + arcPType* RuleDantzig( void ); + +/**< This method returns an arc which violates the dual conditions. It searchs + the arc with most violation of dual conditions in the entire set of real + arcs. It can be used only by the Primal Simplex in the case of networks + with linear costs. */ + +/*--------------------------------------------------------------------------*/ + + arcPType* PRuleFirstEligibleArc( void ); + +/**< This method returns the first found arc which violates the dual conditions + in the case of Primal Simplex, the primal condition in the case of Dual + Simplex. It can be used only in the case of networks with linear costs. */ + +/*--------------------------------------------------------------------------*/ + + arcDType* DRuleFirstEligibleArc( void ); + +/**< This method returns the first found arc which violates the dual conditions + in the case of Primal Simplex, the primal condition in the case of Dual + Simplex. It can be used only in the case of networks with linear costs. */ + +/*--------------------------------------------------------------------------*/ + + arcPType* RulePrimalCandidateListPivot( void ); + +/**< This method returns an arc which violates the dual conditions. It searches + the arc with most violation of dual conditions in a small set of candidate + arcs. In every iteration the method rebuilds this set of arcs executing three + phases: + + - in the first phase it analyzes the remaining arcs and delete the arcs which + don't violate the dual condition any more; + + - in the second phase it tries to fill the set, so it searchs other arcs + which violate the dual condition: the set of arcs is divided into "buckets" + which are searched sequentially until the candidate list is full; the + last visited bucket is retained, and the search is restarted from that + one at later iterations + + - in the third phase the small set of candidate arcs is ordered according + to the violation of dual condition by the method SortPrimalCandidateList() + using an implementation of the algorithm "quicksort". + + At last the method returns the first arc in the ordered small set. If the + arc doesn't exist (the set is empty), it returns NULL. */ + +/*--------------------------------------------------------------------------*/ + + inline void InitializePrimalCandidateList( void ); + +/**< Method to initialize some important details for Primal Candidate List + Rule. */ + +/*--------------------------------------------------------------------------*/ + + inline void SortPrimalCandidateList( Index min , Index max ); + +/**< Method to order the little set of candidate arcs according to + infeasibility of dual conditions. It implements the "quicksort" + algorithm. */ + +/*--------------------------------------------------------------------------*/ + + arcDType* RuleDualCandidateListPivot( void ); + +/**< Similar to RulePrimalCandidateListPivot() for the Dual Simplex. */ + +/*--------------------------------------------------------------------------*/ + + inline void InitializeDualCandidateList( void ); + +/**< Method to initialize some important details for Dual Candidate List Rule. + */ + +/*--------------------------------------------------------------------------*/ + + inline void SortDualCandidateList( Index min , Index max ); + +/**< Similar to SortPrimalCandidateList() for the Dual Simplex. */ + +/*--------------------------------------------------------------------------*/ + + template<class N, class RCT> + inline void AddPotential( N *r , RCT delta ); + +/**< Method to quickly update the dual solutions. During the change of the + base, the potential of node "k2" (deepest node in T of incoming arc "k", + and new root of T2) changes according to the new structure of T. In fact, + the precedessor of "k2" changes: now the predecessor of "k2"is "k1" (the + other node of the incoming arc "k"). This change of predecessor causes a + change of potential of "k2". The change of potential of "k2" causes the + changes of potential of all nodes of T2. This method computes the change + of potential of "k2" and applies it on all the nodes of T2. */ + +/*--------------------------------------------------------------------------*/ + + template<class N> + inline void ComputePotential( N *r ); + +/**< Method to update the dual solutions. It computes all the potential of the + nodes of the subtree which has r as root. */ + +/*--------------------------------------------------------------------------*/ + + inline void ResetWhenInT2( void ); + +/**< Method to order the small set of candidate arcs according to dual + infeasibility. It implements the algorithm "quicksort". */ + +/*--------------------------------------------------------------------------*/ + + void CreateInitialPModifiedBalanceVector( void ); + +/**< Method to initialize the vector of modified balance for the postvisit on + T in the Primal Simplex data structure. */ + +/*--------------------------------------------------------------------------*/ + + void PostPVisit( nodePType *r ); + +/**< Method to calculate the flow on the basic arcs with the Primal Simplex's + data structure. It uses the set of the upper bound arcs, the construction + of a modified balance vector and the postvisit on T. */ + +/*--------------------------------------------------------------------------*/ + + void BalanceFlow( nodePType *r ); + +/**< This method works after the method PostPVisit (in a Primal context). It + restores primal admissimibility on the r's subtree. It starts from the leaf + of the subtree and goes up to the root, using the method AdjustFlow. */ + +/*--------------------------------------------------------------------------*/ + + void AdjustFlow( nodePType *r ); + +/**< This method checks the primal admissibility of the r's basic entering arc. + If it is out of bounds, the method removes it from the tree (and keeps the + relative dummy arc) and push flow in the cycle (some tree's arc and the old + entering arc) to restores the right balances of the node. */ + +/*--------------------------------------------------------------------------*/ + + void CreateInitialDModifiedBalanceVector( void ); + +/**< Method to initialize the vector of modified balance for the postvisit on + T in the Dual Simplex data structure. */ + +/*--------------------------------------------------------------------------*/ + + void PostDVisit( nodeDType *r ); + +/**< Method to calculate the flow on the basic arcs with the Dual Simplex's data + structure, using the set of the upper bound arcs, the construction of a + modified balance vector and the postvisit on T. */ + +/*--------------------------------------------------------------------------*/ + + template<class N, class A> + inline N* Father( N *n, A *a ); + +/**< Method to find the predecessor of the node in the tree. */ + +/*--------------------------------------------------------------------------*/ + + #if(QUADRATICCOST) + template<class A> + inline FONumber ReductCost( A *a ); + #else + template<class A> + inline CNumber ReductCost( A *a ); + #endif + +/**< Method to calculate the reduct cost of the arc. */ + +/*--------------------------------------------------------------------------*/ + + inline FONumber GetFO( void ); + +/**< Method to calculate the temporary (or the final) objective function + value. */ + +/*--------------------------------------------------------------------------*/ + + void PrintPNode( nodePType *nodo ); + +/**< Method to print the "name" of the node in the Primal Simplex. */ + +/*--------------------------------------------------------------------------*/ + + void PrintPArc( arcPType *arc ); + +/**< Method to print the "name" of the arc in the Primal Simplex. */ + +/*--------------------------------------------------------------------------*/ + + void PrintDNode( nodeDType *nodo ); + +/**< Method to print the "name" of the node in the Dual Simplex. */ + +/*--------------------------------------------------------------------------*/ + + void PrintDArc( arcDType *arc ); + +/**< Method to print the "name" of the arc in the Dual Simplex. */ + +/*--------------------------------------------------------------------------*/ + + nodePType* RecoverPNode( Index ind ); + +/**< Method to find a node (in the Primal Simplex) using its index. */ + +/*--------------------------------------------------------------------------*/ + + arcPType* RecoverPArc( nodePType *tail , nodePType *head ); + +/**< Method to find an arc (in the Primal Simplex) using 2 pointers to tail + node and head node. */ + +/*--------------------------------------------------------------------------*/ + + nodeDType* RecoverDNode( Index ind ); + +/**< Method to find a node (in the Dual Simplex) using its index. */ + +/*--------------------------------------------------------------------------*/ + + arcDType* RecoverDArc( nodeDType *tail , nodeDType *head ); + +/**< Method to find an arc (in the Dual Simplex) using 2 pointers to tail + node and head node. */ + +/*--------------------------------------------------------------------------*/ + + void infoPNode( nodePType *node , int tab ); + +/**< Method to print some information of the node (in the Primal Simplex). */ + +/*--------------------------------------------------------------------------*/ + + void infoPArc( arcPType *arc , int ind , int tab ); + +/**< Method to print some information of the arc (in the Primal Simplex). */ + +/*--------------------------------------------------------------------------*/ + + void infoDNode( nodeDType *node , int tab ); + +/**< Method to print some information of the node (in the Dual Simplex). */ + +/*--------------------------------------------------------------------------*/ + + void infoDArc( arcDType *arc , int ind , int tab ); + +/**< Method to print some information of the arc (in the Dual Simplex). */ + +/*--------------------------------------------------------------------------*/ + + void ShowSituation( int tab ); + +/**< Method to show the actual complete situation. */ + +/*--------------------------------------------------------------------------*/ + + }; // end( class MCFSimplex ) + +/* @} end( group( MCFSimplex_CLASSES ) ) */ + +#endif /* MCFSimplex.h included */ + +/*-------------------------------------------------------------------------*/ +/*-------------------inline methods implementation-------------------------*/ +/*-------------------------------------------------------------------------*/ + +inline MCFClass::Index MCFSimplex::MCFSNde( MCFClass::cIndex i ) +{ + if( usePrimalSimplex ) + return( Index( ( (arcsP + i)->tail - nodesP + 1 ) - USENAME0 ) ); + else + return( Index( ( (arcsD + i)->tail - nodesD + 1 ) - USENAME0 ) ); + } + +/*-------------------------------------------------------------------------*/ + +inline MCFClass::Index MCFSimplex::MCFENde( MCFClass::cIndex i ) +{ + if( usePrimalSimplex ) + return( Index( ( (arcsP + i)->head - nodesP + 1 ) - USENAME0 ) ); + else + return( Index( ( (arcsD + i)->head - nodesD + 1 ) - USENAME0 ) ); + } + +/*-------------------------------------------------------------------------*/ + +inline MCFClass::CNumber MCFSimplex::MCFCost( MCFClass::cIndex i ) +{ + if( usePrimalSimplex ) + return( (arcsP + i)->cost ); + else + return( (arcsD + i)->cost ); + } + +/*-------------------------------------------------------------------------*/ + +inline MCFClass::CNumber MCFSimplex::MCFQCoef( MCFClass::cIndex i ) +{ + #if( QUADRATICCOST ) + if( usePrimalSimplex ) + return( (arcsP + i)->quadraticCost ); + else + return( (arcsD + i)->quadraticCost ); + #else + return( 0 ); + #endif + } + +/*-------------------------------------------------------------------------*/ + +inline MCFClass::FNumber MCFSimplex::MCFUCap( MCFClass::cIndex i ) +{ + if( usePrimalSimplex ) + return( (arcsP + i)->upper ); + else + return( (arcsD + i)->upper ); + } + +/*-------------------------------------------------------------------------*/ + +inline MCFClass::FNumber MCFSimplex::MCFDfct( MCFClass::cIndex i ) +{ + if( usePrimalSimplex ) + return( (nodesP + i)->balance ); + else + return( (nodesD + i)->balance ); + } + +/*-------------------------------------------------------------------------*/ + +#if( QUADRATICCOST ) + +template <class A> +inline MCFSimplex::FONumber MCFSimplex::ReductCost( A *a ) +{ + FONumber redc = (a->tail)->potential - (a->head)->potential; + redc = redc + a->cost; + redc = redc + a->quadraticCost * a->flow; + return( redc ); + } + +#else + +template <class A> +inline MCFSimplex::CNumber MCFSimplex::ReductCost( A *a ) +{ + CNumber redc = (a->tail)->potential - (a->head)->potential; + redc = redc + a->cost; + return( redc ); + } + +#endif + +/*-------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) +}; // end( namespace MCFClass_di_unipi_it ) +#endif + +/*-------------------------------------------------------------------------*/ +/*-------------------------------------------------------------------------*/ + +/*-------------------------------------------------------------------------*/ +/*---------------------- End File MCFSimplex.h ----------------------------*/ +/*-------------------------------------------------------------------------*/ diff --git a/impl/MCFSimplex.o b/impl/MCFSimplex.o Binary files differnew file mode 100644 index 0000000..f9378b6 --- /dev/null +++ b/impl/MCFSimplex.o diff --git a/impl/Makefile b/impl/Makefile index ca0b37d..492528d 100644 --- a/impl/Makefile +++ b/impl/Makefile @@ -2,7 +2,8 @@ CC=gcc CPP=g++ BUILD=build PARSER=parser -FLAGS=-Wall -Werror -Wextra -pedantic -lantlr3c -fno-exceptions -lrt +FLAGS= -lantlr3c -fno-exceptions -lrt +#-Wall -Werror -Wextra -pedantic -lantlr3c -fno-exceptions -lrt NORMAL_FLAGS=$(FLAGS) -g -pg OPTIMISED_FLAGS=$(FLAGS) -O3 @@ -16,16 +17,19 @@ check-syntax: $(PARSER) -$(CPP) $(FLAGS) -o /tmp/null -x c++ $(CHK_SOURCES) # the - is to ignore the return code. -$(BUILD)/main: $(BUILD) $(PARSER) main.cpp *.hpp - $(CPP) main.cpp $(PARSER)/*.o -o $(BUILD)/main $(NORMAL_FLAGS) +$(BUILD)/main: $(BUILD) $(PARSER) main.o MCFSimplex.o + $(CPP) main.o MCFSimplex.o $(PARSER)/*.o -o $(BUILD)/main $(NORMAL_FLAGS) -$(BUILD)/fast: $(BUILD) $(PARSER) main.cpp *.hpp - $(CPP) main.cpp $(PARSER)/*.o -o $(BUILD)/fast $(OPTIMISED_FLAGS) +main.o: main.cpp *.hpp *.h + $(CPP) -c main.cpp $(NORMAL_FLAGS) + +MCFSimplex.o: MCFSimplex.cpp *.h + $(CPP) -c MCFSimplex.cpp $(NORMAL_FLAGS) $(PARSER): EquationSystem.g mkdir -p $(PARSER) java -jar ../antlr/antlr-3.4-complete.jar -o $(PARSER) EquationSystem.g - cd $(PARSER); $(CC) *.c -c -lantrl3c + cd $(PARSER); $(CC) *.c -c -lantlr3c $(BUILD): mkdir -p $(BUILD) @@ -37,3 +41,4 @@ test: $(BUILD)/main clean: rm -rf $(BUILD) rm -rf $(PARSER) + rm -rf *.o diff --git a/impl/OPTUtils.h b/impl/OPTUtils.h new file mode 100755 index 0000000..66ad4af --- /dev/null +++ b/impl/OPTUtils.h @@ -0,0 +1,475 @@ +/*--------------------------------------------------------------------------*/ +/*--------------------------- File OPTutils.h ------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @file + * Small classes are provided for: + * - reading the time of a code; + * - generating random numbers. + * + * The classes can be adapted to different environments setting a + * compile-time switch in this file. + * + * Additionally, a function is provided for safely reading parameters + * out of a stream. + * + * \version 1.00 + * + * \date 04 - 10 - 2010 + * + * \author Antonio Frangioni \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * Copyright © 1994 - 2010 by Antonio Frangioni + */ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef __OPTutils + #define __OPTutils /* self-identification: #endif at the end of the file */ + +/*--------------------------------------------------------------------------*/ +/*----------------------------- MACROS -------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup OPTUTILS_MACROS Compile-time switches in OPTutils.h + These macros control how the classes OPTTimers and OPTrand are + implemented; choose the appropriate value for your environment, + or program a new version if no value suits you. + Also, namespaces can be eliminated if they create problems. + @{ */ + + +/*----------------------- OPT_USE_NAMESPACES -------------------------------*/ + +#define OPT_USE_NAMESPACES 0 + +/**< Setting OPT_USE_NAMESPACES == 0 should instruct all codes that use + OPTutils stuff to avoid using namespaces; to start with, the common + namespace OPTutils_di_unipi_it, that contains all the types defined + herein, is *not* defined. */ + +/*---------------------------- OPT_TIMERS ----------------------------------*/ + +#define OPT_TIMERS 5 + +/**< The class OPTtimers is defined below to give an abstract interface to the + different timing routines that are used in different platforms. This is + needed since time-related functions are one of the less standard parts of + the C[++] library. The value of the OPT_TIMERS constant selects among the + different timing routines: + + - 1 = Use the Unix times() routine in sys/times.h + + - 2 = As 1 but uses sys/timeb.h (typical for Microsoft(TM) compilers) + + - 3 = Still use times() of sys/times.h, but return wallclock time + rather than CPU time + + - 4 = As 3 but uses sys/timeb.h (typical for Microsoft(TM) compilers) + + - 5 = return the user time obtained with ANSI C clock() function; this + may result in more accurate running times w.r.t. but may be limited + to ~ 72 hours on systems where ints are 32bits. + + - 6 = Use the Unix gettimeofday() routine of sys/time.h. + + Any unsupported value would simply make the class to report constant + zero as the time. + + The values 1 .. 4 rely on the constant CLK_TCK for converting between + clock ticks and seconds; for the case where the constant is not defined by + the compiler -- should not happen, but it does -- or it is defined in a + wrong way, the constant is re-defined below. */ + +/*---------------------------- OPT_RANDOM ---------------------------------*/ + +#define OPT_RANDOM 1 + +/**< The class OPTrand is defined below to give an abstract interface to the + different random generators that are used in different platforms. This is + needed since random generators are one of the less standard parts of the + C[++] library. The value of the OPT_RANDOM constant selects among the + different timing routines: + + - 0 = an hand-made implementation of a rather good random number generator + is used; note that this assumes that long ints >= 32 bits + + - 1 = standard rand() / srand() pair, common to all C libreries but not + very sophisticated + + - 2 = drand48() / srand48(), common on Unix architectures and pretty good. + + Any unsupported value would simply make the functions to report constant + zero, which is not nice but useful to quickly fix problems if you don't + use random numbers at all. */ + +/*@} -----------------------------------------------------------------------*/ +/*------------------------------ INCLUDES ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( OPT_RANDOM ) + #include <stdlib.h> + + /* For random routines, see OPTrand() below. */ +#endif + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ + +#if( OPT_TIMERS <= 4 ) + #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) ) + #include <sys/times.h> + #else + #include <sys/timeb.h> + #endif + #include <limits.h> +#elif( OPT_TIMERS == 5 ) + #include <time.h> +#elif( OPT_TIMERS == 6 ) + #include <sys/time.h> +#endif + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ + +#include <iostream> +#include <sstream> + +/* For istream and the >> operator, used in DfltdSfInpt(). */ + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ + +/*--------------------------------------------------------------------------*/ +/*--------------------------- NAMESPACE ------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) +namespace OPTtypes_di_unipi_it +{ + /** @namespace OPTtypes_di_unipi_it + The namespace OPTtypes_di_unipi_it is defined to hold all the data + types, constants, classes and functions defined here. It also + comprises the namespace std. */ +#endif + + using namespace std; // I know it's not elegant, but ... + +/*--------------------------------------------------------------------------*/ +/*----------------------------- NULL ---------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef NULL + #define NULL 0 +#endif + +/* Curiously enough, some compilers do not always define NULL; for instance, + sometimes it is defined in stdlib.h. */ + +/*--------------------------------------------------------------------------*/ +/*---------------------------- CLK_TCK -------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef CLK_TCK + #define CLK_TCK 100 +#endif + +/* CLK_TCK is the constant used to transform (integer) time values in seconds + for times()-based routines. Its normal value is 100. */ + +/*--------------------------------------------------------------------------*/ +/*--------------------------- OPT_TIMERS -----------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup OPTtypes_CLASSES Classes in OPTutils.h + @{ */ + +#if( OPT_TIMERS ) + +/** Provides a common interface to the different timing routines that are + available in different platforms. */ + +class OPTtimers { + + public: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + /// constructor of the class + OPTtimers( void ) + { + ReSet(); + } + + /// start the timer + void Start( void ) + { + if( ! ticking ) { + #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 2 ) ) + times( &buff ); + t_u = buff.tms_utime; + t_s = buff.tms_stime; + #elif( ( OPT_TIMERS == 3 ) || ( OPT_TIMERS == 4 ) ) + t_u = times( &buff ); + #elif( OPT_TIMERS == 5 ) + t_u = clock(); + #elif( OPT_TIMERS == 6 ) + struct timeval t; + gettimeofday( &t , NULL ); + t_u = double( t.tv_sec + t.tv_usec * 1e-6 ); + #endif + + ticking = 1; + } + } + + /// stop the timer + void Stop( void ) + { + if( ticking ) { + Read( u , s ); + ticking = 0; + } + } + + /** Return the elapsed time. If the clock is ticking, return the *total* + time since the last Start() without stopping the clock; otherwise, + return the total elapsed time of all the past runs of the clock since + the last ReSet() [see below]. */ + + double Read( void ) + { + double tu = 0; + double ts = 0; + Read( tu , ts ); + return( tu + ts ); + } + + /// As Read( void ) but *adds* user and system time to tu and ts. + void Read( double &tu , double &ts ) + { + if( ticking ) { + #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 2 ) ) + times( &buff ); + tu += ( double( buff.tms_utime - t_u ) ) / double( CLK_TCK ); + ts += ( double( buff.tms_stime - t_s ) ) / double( CLK_TCK ); + #elif( ( OPT_TIMERS == 3 ) || ( OPT_TIMERS == 4 ) ) + tu += ( double( times( &buff ) - t_u ) ) / double( CLK_TCK ); + #elif( OPT_TIMERS == 5 ) + tu += ( double( clock() - t_u ) ) / double( CLOCKS_PER_SEC ); + #elif( OPT_TIMERS == 6 ) + struct timeval t; + gettimeofday( &t , NULL ); + tu += double( t.tv_sec + t.tv_usec * 1e-6 ) - t_u; + #endif + } + else { tu += u; ts += s; } + } + + /// reset the timer + void ReSet( void ) + { + u = s = 0; ticking = 0; + } + + private: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + double u; // elapsed *user* time, in seconds + double s; // elapsed *system* time, in seconds + char ticking; // 1 if the clock is ticking + + #if( ( OPT_TIMERS > 0 ) && ( OPT_TIMERS <= 5 ) ) + clock_t t_u; + + #if( OPT_TIMERS <= 4 ) + struct tms buff; + + #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) ) + clock_t t_s; + #endif + #endif + #elif( OPT_TIMERS == 6 ) + double t_u; + #endif + + }; // end( class OPTtimers ); + +#endif + +/*--------------------------------------------------------------------------*/ +/*------------------------------ OPTrand() ---------------------------------*/ +/*--------------------------------------------------------------------------*/ + +/** Provide a common interface to the different random generators that are + available in different platforms. */ + +class OPTrand { + + public: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + /// constructor of the class + OPTrand( void ) + { + #if( OPT_RANDOM == 0 ) + A[ 0 ] = -1; + srand( long( 123456789 ) ); + #else + OPTrand::srand( long( 1 ) ); + #endif + } + + /** Returns a random number uniformly distributed in [0, 1). + \note each object of class OPTrand has its own sequence, so that + multiple OPTrand objects being used within the same program do not + interfere with each other (as opposed to what C random routines + would do). */ + + double rand( void ) + { + #if( OPT_RANDOM == 0 ) + long nmbr = *(gb_fptr--); + if( nmbr < 0 ) + nmbr = gb_flip_cycle(); + return( double( nmbr ) / double( (unsigned long)0x80000000 ) ); + #elif( OPT_RANDOM == 1 ) + ::srand( myseed ); + myseed = ::rand(); + return( double( myseed ) / double( RAND_MAX ) ); + #elif( OPT_RANDOM == 2 ) + return( erand48( myseed ) ); + #else + return( 0 ); // random, eh? + #endif + } + + /// Seeds the random generator for this instance of OPTrand. + void srand( long seed ) + { + #if( OPT_RANDOM == 0 ) + long prev = seed , next = 1; + seed = prev = mod_diff( prev , 0 ); + A[ 55 ] = prev; + + for( long i = 21 ; i ; i = ( i + 21 ) % 55 ) { + A[ i ] = next; + next = mod_diff( prev , next ); + if( seed & 1 ) + seed = 0x40000000 + ( seed >> 1 ); + else + seed >>= 1; + + next = mod_diff( next , seed ); + prev = A[ i ]; + } + + gb_flip_cycle(); + gb_flip_cycle(); + gb_flip_cycle(); + gb_flip_cycle(); + gb_flip_cycle(); + #elif( OPT_RANDOM == 1 ) + myseed = int( seed ); + #elif( OPT_RANDOM == 2 ) + long *sp = (long*)( &myseed ); + *sp = seed; // copy higher 32 bits + myseed[ 2 ] = 0x330E; // initialize lower 16 bits + #endif + } + + private: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + #if( OPT_RANDOM == 0 ) + long A[ 56 ]; + long *gb_fptr; + + long mod_diff( long x , long y ) + { + return( ( ( x ) - ( y ) ) & 0x7fffffff ); + } + + long gb_flip_cycle( void ) + { + long *ii, *jj; + for( ii = &A[ 1 ] , jj = &A[ 32 ] ; jj <= &A[ 55 ] ; ii++ , jj++ ) + *ii = mod_diff( *ii , *jj ); + + for( jj = &A[ 1 ] ; ii <= &A[ 55 ] ; ii++ , jj++ ) + *ii = mod_diff( *ii , *jj ); + + gb_fptr = &A[ 54 ]; + + return A[ 55 ]; + } + #elif( OPT_RANDOM == 1 ) + int myseed; + #elif( OPT_RANDOM == 2 ) + unsigned short int myseed[ 3 ]; + #endif + + }; // end( class( OPTrand ) ) + +/* @} end( group( OPTtypes_CLASSES ) ) */ +/*--------------------------------------------------------------------------*/ +/*--------------------------- DfltdSfInpt() --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup OPTtypes_FUNCTIONS Functions in OPTutils.h + @{ */ + +/** Template function for reading parameters from a istream. The function is + "safe" because it works also if the istream is not given, is not be long + enough or contains erroneous things. + + Given a &istream (possibly NULL), DfltdSfInpt() attempts to read Param out + of it, skipping any line that begins with the comment carachter (defaulted + to '#'), any blank line and any line starting with anything that can not + be interpreted as a `T'. If, for any reason, the read operation fails, + then the parameter is given the default value `Dflt'. Otherwise, all the + rest of the line up to the nearest newline ('\n') carachter is flushed. */ + +template<class T> +inline void DfltdSfInpt( istream *iStrm , T &Param , const T Dflt , + const char cmntc = '#' ) +{ + string comm( 1 , cmntc ); + + // the "! ! stream" trick is there to force the compiler to apply the + // stream -> bool conversion, in case it is too dumb to do it by itself + if( iStrm && ( ! ( ! (*iStrm) ) ) ) { + string buf; + + for( ;; ) { + if( ! ( (*iStrm) >> ws ) ) // skip whitespace + break; + + if( ! (*iStrm).good() ) // check stream is OK + break; + + getline( *iStrm , buf ); + + if( ! buf.empty() ) + if( buf.substr( 0 , 1 ).compare( comm ) ) { + stringstream ss; + ss << buf; + if( ! ( ss >> Param ) ) + break; + + return; + } + } + } + + Param = Dflt; + + } // end( DfltdSfInpt ) + +/* @} end( group( OPTtypes_FUNCTIONS ) ) */ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) + }; // end( namespace OPTtypes_di_unipi_it ) +#endif + +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#endif /* OPTutils.h included */ + +/*--------------------------------------------------------------------------*/ +/*---------------------- End File OPTutils.h -------------------------------*/ +/*--------------------------------------------------------------------------*/ diff --git a/impl/Operator.hpp b/impl/Operator.hpp index 64ef096..4a573b8 100644 --- a/impl/Operator.hpp +++ b/impl/Operator.hpp @@ -135,6 +135,88 @@ struct Guard : public Operator<Domain> { } }; +#include "MCFSimplex.h" + +template<typename Domain> +struct MinCostFlow : public Operator<Domain> { + MinCostFlow(const std::vector<Domain>& supplies, const std::vector<std::pair<int,int> > arcs) + : _supplies(supplies), _arcs(arcs), _solver(0,0) { + MCFClass::FNumber* deficits = new MCFClass::FNumber[supplies.size()]; + MCFClass::Index* starts = new MCFClass::Index[arcs.size()]; + MCFClass::Index* ends = new MCFClass::Index[arcs.size()]; + + for (int i = 0, size = supplies.size(); i < size; ++i) { + deficits[i] = -supplies[i].template as<MCFClass::FNumber>(); + } + for (int i = 0, size = arcs.size(); i < size; ++i) { + starts[i] = arcs[i].first; + ends[i] = arcs[i].second; + } + + _solver.LoadNet(supplies.size(), supplies.size(), // max nodes/arcs + supplies.size(), supplies.size(), // current nodes/arcs + NULL, NULL, // arcs have inf cap, zero cost (to begin) + deficits, // deficits for each node + starts, ends); // start/end of each edge + + delete[] deficits; + delete[] starts; + delete[] ends; + } + Domain eval (const std::vector<Domain>& costs) const { + assert(costs.size() == _arcs.size()); + if (costs.size() == 0) + return 0; + for (int i = 0, size = costs.size(); i < size; ++i) { + _solver.ChgCost(i, costs[i].template as<MCFClass::CNumber>()); + } + _solver.SolveMCF(); + if (_solver.MCFGetStatus() == MCFClass::kUnfeasible){ + return -infinity<Domain>(); + } else if (_solver.MCFGetFO() == MCFClass::Inf<MCFClass::FONumber>()) { + return infinity<Domain>(); + } else if (_solver.MCFGetFO() == -MCFClass::Inf<MCFClass::FONumber>()) { + return -infinity<Domain>(); + } else { + return _solver.MCFGetFO(); + } + } + void print(std::ostream& cout) const { + std::string supplyString = "["; + for (int i = 0, size = _supplies.size(); i < size; ++i) { + if (i > 0) + supplyString += ","; + std::stringstream stream; + stream << _supplies[i]; + supplyString += stream.str(); + } + supplyString += ']'; + + std::string arcString = "["; + for (int i = 0, size = _arcs.size(); i < size; ++i) { + if (i > 0) + arcString += ","; + { + std::stringstream stream; + stream << _arcs[i].first; + arcString += stream.str() + ":"; + } + { + std::stringstream stream; + stream << _arcs[i].second; + arcString += stream.str(); + } + } + arcString += ']'; + + cout << "MCF<" << supplyString << "," << arcString << ">"; + } +private: + std::vector<Domain> _supplies; + std::vector<std::pair<int,int> > _arcs; + mutable MCFSimplex _solver; +}; + /*#include "TemplateConstraintMatrix.hpp" template<typename Domain> diff --git a/impl/main.cpp b/impl/main.cpp index be9cc82..0daa72a 100644 --- a/impl/main.cpp +++ b/impl/main.cpp @@ -20,6 +20,59 @@ extern "C" { using namespace std; template<typename T> +vector<T> toSupplies(pANTLR3_BASE_TREE node) { + ANTLR3_UINT32 num = node->getChildCount(node); + T output; + vector<T> supplies; + pANTLR3_BASE_TREE childNode; + + for (int i = 0; i < num; ++i) { + childNode = (pANTLR3_BASE_TREE) node->getChild(node, i); + stringstream stream((char*)childNode->getText(childNode)->chars); + assert(stream >> output); + supplies.push_back(output); + } + + return supplies; +} + +pair<int,int> toArc(pANTLR3_BASE_TREE node) { + pair<int,int> result; + int output; + pANTLR3_BASE_TREE childNode; + + ANTLR3_UINT32 num = node->getChildCount(node); + assert(num == 2); + + { + childNode = (pANTLR3_BASE_TREE) node->getChild(node, 0); + stringstream stream((char*)childNode->getText(childNode)->chars); + assert(stream >> output); + result.first = output; + } + + { + childNode = (pANTLR3_BASE_TREE) node->getChild(node, 1); + stringstream stream((char*)childNode->getText(childNode)->chars); + assert(stream >> output); + result.second = output; + } + + return result; +} + +vector<pair<int,int> > toArcs(pANTLR3_BASE_TREE node) { + vector<pair<int,int> > arcs; + ANTLR3_UINT32 num = node->getChildCount(node); + pANTLR3_BASE_TREE childNode; + for (int i = 0; i < num; ++i) { + childNode = (pANTLR3_BASE_TREE) node->getChild(node, i); + arcs.push_back(toArc(childNode)); + } + return arcs; +} + +template<typename T> Expression<T>& treeToExpression(pANTLR3_BASE_TREE node, EquationSystem<T>& system) { ANTLR3_UINT32 num = node->getChildCount(node); string name = (char*) node->getText(node)->chars; @@ -39,7 +92,23 @@ Expression<T>& treeToExpression(pANTLR3_BASE_TREE node, EquationSystem<T>& syste } } - // anything that's not a constant/variable + if (name == "MCF") { + pANTLR3_BASE_TREE childNode; + + childNode = (pANTLR3_BASE_TREE) node->getChild(node, 0); + vector<T> supplies = toSupplies<T>(childNode); + childNode = (pANTLR3_BASE_TREE) node->getChild(node, 1); + vector<pair<int,int> > arcs = toArcs(childNode); + + vector<Expression<T>*> args; + for (unsigned int i = 2; i < num; ++i) { + childNode = (pANTLR3_BASE_TREE) node->getChild(node, i); + args.push_back(&treeToExpression(childNode, system)); + } + return system.expression(new MinCostFlow<T>(supplies, arcs), args); + } + + // anything that's not a constant/variable, or an MCF expr vector<Expression<T>*> args; pANTLR3_BASE_TREE childNode; for (unsigned int i = 0; i < num; ++i) { @@ -70,7 +139,7 @@ Expression<T>& treeToExpression(pANTLR3_BASE_TREE node, EquationSystem<T>& syste } else if (name == "GUARD" || name == "guard") { op = new Guard<T>(); } else { - std::cout << "throw exception" << *(char*)NULL; + assert(false); //throw "Parse error: Unknown operator"; } return system.expression(op, args); |