/*** 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