/*** Copyright 1989 Norbert Schlenker.  All rights reserved.

  *** This software is distributed subject to the following provisions:

  ***    - this notice may not be removed;

  ***    - you may modify the source code, as long as redistributed

  ***      versions have their modifications clearly marked;

  ***    - no charge, other than a nominal copying fee, may be made

  ***      when providing copies of this source code to others;

  ***    - if this source code is used as part of a product which is

  ***      distributed only as a binary, a copy of this source code

  ***      must be included in the distribution.

  ***

  *** Unlike the GNU GPL, use of this code does not obligate you to

  *** disclose your own proprietary source code.

  *** The author of this software provides no warranty, express or

  *** implied, as to its utility or correctness.  That said, reports

  *** of bugs or compatibility problems will be gladly received by

  *** nfs@princeton.edu, and fixes will be attempted.

  ***/

 /*** netgen - C version of the standard NETGEN network generator 

  ***          This program is a functional equivalent of the

  ***          standard network generator NETGEN described in:

  ***		Klingman, D., A. Napier, and J. Stutz, "NETGEN:  A Program

  ***		  for Generating Large Scale Capacitated Assignment,

  ***		  Transportation, and Minimum Cost Flow Network Problems",

  ***		  Management Science 20, 5, 814-821 (1974)

  ***

  ***	      This software provides a number of interfaces for use by

  ***	      network solvers.  Standard call interfaces are supplied for

  ***	      use by (Unix) C and Fortran solvers, with generation parameters

  ***	      passed into the generator and the flow network passed back to

  ***	      the solver via large external (COMMON in Fortran) arrays.

  ***	      For the DIMACS challenge, this code will produce output files

  ***	      in the appropriate format for later reading by solvers.

  ***          Undefine the symbol DIMACS when using the call interface.

  ***

  ***          The generator produces exact duplicates of the networks

  ***          made by the Fortran code (even though that means bugs

  ***          are being perpetuated). It is faster by orders of magnitude

  ***          in generating large networks, primarily by imposing the

  ***          notion of the abstract data type INDEX_LIST and implementing

  ***          the type operations in a reasonably efficient fashion.

  ***/

 /*** Generates transportation problems if:

  ***	SOURCES+SINKS == NODES && TSOURCES == TSINKS == 0

  ***

  *** Generates assignment problems if above conditions are satisfied, and:

  ***	SOURCES == SINKS && SUPPLY == SOURCES

  ***

  *** Generates maximum flow problems if not an assignment problem and:

  ***	MINCOST == MAXCOST == 1

  *** Implementation notes:

  ***

  ***	This file contains both a Fortran and a C interface. The

  ***	Fortran interface is suffixed with an underscore to make it

  ***	callable in the normal fashion from Fortran (a Unix convention).

  ***

  ***    Because Fortran has no facility for pointers, the common arrays

  ***    are statically allocated.  Static allocation has nothing to recommend

  ***    it except for the need for a Fortran interface.

  ***

  ***    This software expects input parameters to be long integers

  ***    (in the sense of C); that means no INTEGER*2 from Fortran callers.

  ***

  ***	Compiling with -DDIMACS produces a program that reads problem

  ***	parameters, generates the appropriate problem, and prints it.

  ***

  ***	Compiling with -DDEBUG produces code with externally visible

  ***	procedure names, useful for debugging and profiling.

  ***/

 /*** System interfaces */

 #include <stdio.h>

 /*** Public interfaces */

 #define ALLOCATE_NETWORK

 #include "netgen.h"

 #define PROBLEM_PARMS 13		/* aliases for generation parameters */

 #define NODES	    parms[0]		/* number of nodes */

 #define SOURCES     parms[1]		/* number of sources (including transshipment) */

 #define SINKS	    parms[2]		/* number of sinks (including transshipment) */

 #define DENSITY     parms[3]		/* number of (requested) arcs */

 #define MINCOST     parms[4]		/* minimum cost of arcs */

 #define MAXCOST     parms[5]		/* maximum cost of arcs */

 #define SUPPLY	    parms[6]		/* total supply */

 #define TSOURCES    parms[7]		/* transshipment sources */

 #define TSINKS	    parms[8]		/* transshipment sinks */

 #define HICOST	    parms[9]		/* percent of skeleton arcs given maximum cost */

 #define CAPACITATED parms[10]		/* percent of arcs to be capacitated */

 #define MINCAP	    parms[11]		/* minimum capacity for capacitated arcs */

 #define MAXCAP	    parms[12]		/* maximum capacity for capacitated arcs */

 /*** Private interfaces */

 #ifdef DEBUG

 #define PRIVATE

 #else

 #define PRIVATE static

 #endif

 #ifdef __STDC__

 PRIVATE void create_supply(NODE, CAPACITY); /* create supply nodes */

 PRIVATE void create_assignment(long*);	/* create assignment problem */

 PRIVATE void sort_skeleton(int);	/* sorts skeleton chains */

 PRIVATE void pick_head(long*, int, NODE); /* choose destination nodes for rubbish arcs */

 PRIVATE void error_exit(long);		/* print error message and exit */

 #else

 PRIVATE void create_supply();		/* create supply nodes */

 PRIVATE void create_assignment();	/* create assignment problem */

 PRIVATE void sort_skeleton();		/* sorts skeleton chains */

 PRIVATE void pick_head();		/* chooses destination nodes for rubbish arcs */

 PRIVATE void error_exit();		/* print error message and exit */

 #endif

 /*** Private variables */

 static NODE nodes_left;

 static ARC arc_count;

 static NODE pred[MAXARCS];

 static NODE head[MAXARCS];

 static NODE tail[MAXARCS];

 /*** Local macros */

 #define MIN(x, y) ((x) < (y) ? (x) : (y))

 #define MAX(x, y) ((x) > (y) ? (x) : (y))

 #define SAVE_ARC(tail, head, cost, capacity)	/* records an arc where our caller can get it */ \

   {				\

     FROM[arc_count] = tail;	\

     TO  [arc_count] = head;	\

     C   [arc_count] = cost;	\

     U   [arc_count] = capacity; \

     arc_count++;		\

   }

 /*** Fortran callable interface routine */

 void netgen_(seed, parms, generated_nodes, generated_arcs)

 long* seed;			/* pointer to random seed */

 long parms[PROBLEM_PARMS];	/* problem parameters */

 long* generated_nodes;		/* number of generated nodes */

 long* generated_arcs;		/* number of generated arcs */

 {

   *generated_nodes = NODES;

   if ((*generated_arcs = netgen(*seed, parms)) < 0)

     error_exit(*generated_arcs);

 }

 /*** C callable interface routine */

 ARC netgen(seed, parms)

 long seed;			/* random seed */

 long parms[];			/* problem parameters */

 {

   register NODE i,j,k;

   NODE source;

   NODE node;

   NODE sinks_per_source;

   NODE* sinks;

   NODE it;

   int chain_length;

   COST cost;

   CAPACITY cap;

   INDEX_LIST handle;

   int supply_per_sink;

   int partial_supply;

   int sort_count;

 /*** Perform sanity checks on the input */

   if (seed <= 0)

     return BAD_SEED;

   if (NODES > MAXNODES || DENSITY > MAXARCS)

     return TOO_BIG;

   if ((NODES <= 0) ||

       (NODES > DENSITY) ||

       (SOURCES <= 0) ||

       (SINKS <= 0) ||

       (SOURCES + SINKS > NODES) ||

       (MINCOST > MAXCOST) ||

       (SUPPLY < SOURCES) ||

       (TSOURCES > SOURCES) ||

       (TSINKS > SINKS) ||

       (HICOST < 0 || HICOST > 100) ||

       (CAPACITATED < 0 || CAPACITATED > 100) ||

       (MINCAP > MAXCAP))

     return BAD_PARMS;

 /*** Do a little bit of setting up. */

   set_random(seed);

   arc_count = 0;

   nodes_left = NODES - SINKS + TSINKS;

   if ((SOURCES-TSOURCES)+(SINKS-TSINKS) == NODES &&

       (SOURCES-TSOURCES) == (SINKS-TSINKS) &&

        SOURCES == SUPPLY) {

     create_assignment(parms);

     return arc_count;

   }

   (void)memset((void *)B, 0, sizeof(B));/* set supplies and demands to zero */

   create_supply((NODE)SOURCES, (CAPACITY)SUPPLY);

 /*** Form most of the network skeleton.  First, 60% of the transshipment

  *** nodes are divided evenly among the various sources; the remainder

  *** are chained onto the end of the chains belonging to random sources.

  ***/

   for (i = 1; i <= SOURCES; i++)	/* point SOURCES at themselves */

     pred[i] = i;

   handle = make_index_list((INDEX)(SOURCES + 1), (INDEX)(NODES - SINKS));

   source = 1;

   for (i = NODES-SOURCES-SINKS; i > (4*(NODES-SOURCES-SINKS)+9)/10; i--) {

     node = choose_index(handle, (INDEX)ng_random(1L, (long)index_size(handle)));

     pred[node] = pred[source];

     pred[source] = node;

     if (++source > SOURCES)

       source = 1;

   }

   for ( ; i > 0; --i) {

     node = choose_index(handle, (INDEX)ng_random(1L, (long)index_size(handle)));

     source = ng_random(1L, SOURCES);

     pred[node] = pred[source];

     pred[source] = node;

   }

   free_index_list(handle);

 /*** For each source chain, hook it to an "appropriate" number of sinks,

  *** place capacities and costs on the skeleton edges, and then call

  *** pick_head to add a bunch of rubbish edges at each node on the chain.

  ***/

   for (source = 1; source <= SOURCES; source++) {

     sort_count = 0;

     node = pred[source];

     while (node != source) {

       sort_count++;

       head[sort_count] = node;

       node = tail[sort_count] = pred[node];

     }

     if ((NODES-SOURCES-SINKS) == 0)

       sinks_per_source = SINKS/SOURCES + 1;

     else

 /* changing to handle overflows with large n; Mar 18 -- jc */

       sinks_per_source = ((double) 2*sort_count*SINKS) / ((double) NODES-SOURCES-SINKS);

     sinks_per_source = MAX(2, MIN(sinks_per_source, SINKS));

     sinks = (NODE*) malloc(sinks_per_source * sizeof(NODE));

     handle = make_index_list((INDEX)(NODES - SINKS), (INDEX)(NODES - 1));

     for (i = 0; i < sinks_per_source; i++) {

       sinks[i] = choose_index(handle, (INDEX)ng_random(1L, (long)index_size(handle)));

     }

     if (source == SOURCES && index_size(handle) > 0) {

       sinks = (NODE*) realloc((void *)sinks, (sinks_per_source + index_size(handle)) * sizeof(NODE));

       while (index_size(handle) > 0) {

 	j = choose_index(handle, 1);

 	if (B[j] == 0)

 	  sinks[sinks_per_source++] = j;

       }

     }

     free_index_list(handle);

     chain_length = sort_count;

     supply_per_sink = B[source-1] / sinks_per_source;

     k = pred[source];

     for (i = 0; i < sinks_per_source; i++) {

       sort_count++;

       partial_supply = ng_random(1L, (long)supply_per_sink);

       j = ng_random(0L, (long)sinks_per_source - 1);

       tail[sort_count] = k;

       head[sort_count] = sinks[i] + 1;

       B[sinks[i]] -= partial_supply;

       B[sinks[j]] -= (supply_per_sink - partial_supply);

       k = source;

       for (j = ng_random(1L, (long)chain_length); j > 0; j--)

 	k = pred[k];

     }

     B[sinks[0]] -= (B[source-1] % sinks_per_source);

     free((void *)sinks);

     sort_skeleton(sort_count);

     tail[sort_count+1] = 0;

     for (i = 1; i <= sort_count; ) {

       handle = make_index_list((INDEX)(SOURCES - TSOURCES + 1), (INDEX)NODES);

       remove_index(handle, (INDEX)tail[i]);

       it = tail[i];

       while (it == tail[i]) {

 	remove_index(handle, (INDEX)head[i]);

 	cap = SUPPLY;

 	if (ng_random(1L, 100L) <= CAPACITATED)

 	  cap = MAX(B[source-1], MINCAP);

 	cost = MAXCOST;

 	if (ng_random(1L, 100L) > HICOST)

 	  cost = ng_random(MINCOST, MAXCOST);

 	SAVE_ARC(it,head[i],cost,cap);

 	i++;

       }

       pick_head(parms, handle, it);

       free_index_list(handle);

     }

   }

 /*** Add more rubbish edges out of the transshipment sinks. */

   for (i = NODES - SINKS + 1; i <= NODES - SINKS + TSINKS; i++) {

     handle = make_index_list((INDEX)(SOURCES - TSOURCES + 1), (INDEX)NODES);

     remove_index(handle, (INDEX)i);

     pick_head(parms, handle, i);

     free_index_list(handle);

   }

   return arc_count;

 }

 PRIVATE void create_supply(sources, supply)

 NODE sources;

 CAPACITY supply;

 {

   CAPACITY supply_per_source = supply / sources;

   CAPACITY partial_supply;

   NODE i;

   for (i = 0; i < sources; i++) {

     B[i] += (partial_supply = ng_random(1L, (long)supply_per_source));

     B[ng_random(0L, (long)(sources - 1))] += supply_per_source - partial_supply;

   }

   B[ng_random(0L, (long)(sources - 1))] += supply % sources;

 }

 PRIVATE void create_assignment(parms)

 long parms[];

 {

   INDEX_LIST skeleton, handle;

   INDEX index;

   NODE source;

   for (source = 0; source < NODES/2; source++)

     B[source] = 1;

   for ( ; source < NODES; source++)

     B[source] = -1;

   skeleton = make_index_list((INDEX)(SOURCES + 1), (INDEX)NODES);

   for (source = 1; source <= NODES/2; source++) {

     index = choose_index(skeleton, (INDEX)ng_random(1L, (long)index_size(skeleton)));

     SAVE_ARC(source, index, ng_random(MINCOST, MAXCOST), 1);

     handle = make_index_list((INDEX)(SOURCES + 1), (INDEX)NODES);

     remove_index(handle, index);

     pick_head(parms, handle, source);

     free_index_list(handle);

   }

   free_index_list(skeleton);

 }

 PRIVATE void sort_skeleton(sort_count) 		/* Shell sort */

 int sort_count;

 {

   int m,i,j,k;

   int temp;

   m = sort_count;

   while ((m /= 2) != 0) {

     k = sort_count - m;

     for (j = 1; j <= k; j++) {

       i = j;

       while (i >= 1 && tail[i] > tail[i+m]) {

 	temp = tail[i];

 	tail[i] = tail[i+m];

 	tail[i+m] = temp;

 	temp = head[i];

 	head[i] = head[i+m];

 	head[i+m] = temp;

 	i -= m;

       }

     }

   }

 }

 PRIVATE void pick_head(parms, handle, desired_tail)

 long parms[];

 INDEX_LIST handle;

 NODE desired_tail;

 {

   NODE non_sources = NODES - SOURCES + TSOURCES;

 /* changing Aug 29 -- jc

   ARC remaining_arcs = DENSITY - arc_count;

 */

   int  remaining_arcs = (int) DENSITY - (int) arc_count;

   INDEX index;

   int limit;

   long upper_bound;

   CAPACITY cap;

 /* changing Aug 29 -- jc

 */

   nodes_left--;

   if ((2 * (int) nodes_left) >= (int) remaining_arcs)

     return;

   if ((remaining_arcs + non_sources - pseudo_size(handle) - 1) / (nodes_left + 1) >= non_sources - 1) {

     limit = non_sources;

   } else {

     upper_bound = 2 * (remaining_arcs / (nodes_left + 1) - 1);

     do {

       limit = ng_random(1L, upper_bound);

       if (nodes_left == 0)

 	limit = remaining_arcs;

 /* changing to handle overflows with large n; Mar 18 -- jc */

     } while ( ((double) nodes_left * (non_sources - 1)) < ((double) remaining_arcs - limit));

   }

   for ( ; limit > 0; limit--) {

     index = choose_index(handle, (INDEX)ng_random(1L, (long)pseudo_size(handle)));

     cap = SUPPLY;

     if (ng_random(1L, 100L) <= CAPACITATED)

       cap = ng_random(MINCAP, MAXCAP);

 /* adding Aug 29 -- jc */

 if ((1 <= index) && (index <= NODES)) {

     SAVE_ARC(desired_tail, index, ng_random(MINCOST, MAXCOST), cap);

     }

   }

 }

 /*** Print an appropriate error message and then exit with a nonzero code. */

 PRIVATE void error_exit(rc)

 long rc;

 {

   switch (rc) {

     case BAD_SEED:

       fprintf(stderr, "NETGEN requires a positive random seed\n");

       break;

     case TOO_BIG:

       fprintf(stderr, "Problem too large for generator\n");

       break;

     case BAD_PARMS:

       fprintf(stderr, "Inconsistent parameter settings - check the input\n");

       break;

     case ALLOCATION_FAILURE:

       fprintf(stderr, "Memory allocation failure\n");

       break;

     default:

       fprintf(stderr, "Internal error\n");

       break;

   }

   exit(1000 - (int)rc);

 }

 #ifdef DIMACS			/* generates network on standard output */

 #define READ(v) 		     	/* read one variable using scanf */	\

 	switch( scanf("%ld", &v) ) {						\

 		case 1:								\

 			break;							\

 		default:							\

 			exit(0);						\

 		}

 int main()

 {

   long seed;

   long problem;

   long parms[PROBLEM_PARMS];

   long arcs;

   int i;

 /*** Read problem parameters and generate networks */

   while (1) {

     READ(seed);

     if (seed <= 0) exit(0);

     READ(problem);

     if (problem <= 0) exit(0);

     for (i = 0; i < PROBLEM_PARMS; i++)

       READ(parms[i]);

     printf("c NETGEN flow network generator (C version)\n");

     printf("c  Problem %2ld input parameters\n", problem);

     printf("c  ---------------------------\n");

     printf("c   Random seed:          %10ld\n",   seed);

     printf("c   Number of nodes:      %10ld\n",   NODES);

     printf("c   Source nodes:         %10ld\n",   SOURCES);

     printf("c   Sink nodes:           %10ld\n",   SINKS);

     printf("c   Number of arcs:       %10ld\n",   DENSITY);

     printf("c   Minimum arc cost:     %10ld\n",   MINCOST);

     printf("c   Maximum arc cost:     %10ld\n",   MAXCOST);

     printf("c   Total supply:         %10ld\n",   SUPPLY);

     printf("c   Transshipment -\n");

     printf("c     Sources:            %10ld\n",   TSOURCES);

     printf("c     Sinks:              %10ld\n",   TSINKS);

     printf("c   Skeleton arcs -\n");

     printf("c     With max cost:      %10ld%%\n", HICOST);

     printf("c     Capacitated:        %10ld%%\n", CAPACITATED);

     printf("c   Minimum arc capacity: %10ld\n",   MINCAP);

     printf("c   Maximum arc capacity: %10ld\n",   MAXCAP);

     if ((arcs = netgen(seed, parms)) < 0)

       error_exit(arcs);

     if ((SOURCES-TSOURCES)+(SINKS-TSINKS) == NODES &&

 	(SOURCES-TSOURCES) == (SINKS-TSINKS) &&

 	 SOURCES == SUPPLY) {

       printf("c\n");

       printf("c  *** Assignment ***\n");

       printf("c\n");

       printf("p asn %ld %ld\n", NODES, arcs);

       for (i = 0; i < NODES; i++) {

 	if (B[i] > 0)

 	  printf("n %ld\n", i + 1);

       }

       for (i = 0; i < arcs; i++) {

 	printf("a %ld %ld %ld\n", FROM[i], TO[i], C[i]);

       }

     } else

     if (MINCOST == 1 && MAXCOST == 1) {

       printf("c\n");

       printf("c  *** Maximum flow ***\n");

       printf("c\n");

       printf("p max %ld %ld\n", NODES, arcs);

       for (i = 0; i < NODES; i++) {

 	if (B[i] > 0)

 	  printf("n %ld s\n", i + 1);

 	else

 	if (B[i] < 0)

 	  printf("n %ld t\n", i + 1);

       }

       for (i = 0; i < arcs; i++) {

 	printf("a %ld %ld %ld\n", FROM[i], TO[i], U[i]);

       }

     } else {

       printf("c\n");

       printf("c  *** Minimum cost flow ***\n");

       printf("c\n");

       printf("p min %ld %ld\n", NODES, arcs);

       for (i = 0; i < NODES; i++) {

 	if (B[i] != 0)

 	  printf("n %ld %ld\n", i + 1, B[i]);

       }

       for (i = 0; i < arcs; i++) {

 	printf("a %ld %ld %ld %ld %ld\n", FROM[i], TO[i], 0, U[i], C[i]);

       }

     }

   }

   return 0;

 }

 #endif