/*
SFP code: SFP is a Rectilinear Steiner Minimum Tree (RSMT) heuristic.
It accepts inputs in the ISPD'08 format described at
http://www.ispd.cc/contests/08/ispd08rc.html.

Copyright 2019-2022 Texas State University

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Authors: Martin Burtscher, Aarti Kothari, and Alex Fallin

URL: The latest version of this code is available at
https://cs.txstate.edu/~burtscher/research/SFP/.
*/


#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <tuple>
#include <algorithm>
#include <climits>
#include <sys/time.h>


static const int MaxPins = 16384;  // must be a power of 2

using ID = short;  // must be signed (point and edge IDs)
using ctype = int;  // must be signed (coordinates and distances)  // change requires further change below

typedef struct {
  ID src;
  ID dst;
} edge;

static ID source[MaxPins * 2];
static ID destin[MaxPins * 2];
static ctype dist[MaxPins * 2];
static ID adj[MaxPins * 2][8];
static char cnt[MaxPins * 2];


static inline void buildMST(const ID num, const ctype* const __restrict__ x, const ctype* const __restrict__ y, edge* const __restrict__ edges)
{
  // initialize
  ID numItems = num - 1;
  for (ID i = 0; i < numItems; i++) dist[i] = INT_MAX;  // change if ctype changed
  for (ID i = 0; i < numItems; i++) destin[i] = (ID)(i + 1);

  // Prim's MST algorithm
  ID src = 0;
  for (ID cnt = 0; cnt < num - 1; cnt++) {
    int mindj = INT_MAX;

    // update distances
    for (ID j = 0; j < numItems; j++) {
      const ID dst = destin[j];
      const ctype dnew = std::abs(x[src] - x[dst]) + std::abs(y[src] - y[dst]);
      ctype d = dist[j];
      if (d > dnew) {
        d = dnew;
        dist[j] = dnew;
        source[j] = src;
      }
      const int upv = d * (MaxPins * 2) + j;  // tie breaker
      if (mindj > upv) mindj = upv;
    }

    // extract new edge
    const ID j = mindj % (MaxPins * 2);
    src = destin[j];
    edges[cnt].src = source[j];
    edges[cnt].dst = src;

    // remove used element
    numItems--;
    dist[j] = dist[numItems];
    source[j] = source[numItems];
    destin[j] = destin[numItems];
  }
}


static inline bool insertSteinerPoints(ID& num, ctype* const __restrict__ x, ctype* const __restrict__ y, edge* const __restrict__ edges)
{
  const ID top = num;

  // create adjacency lists
  for (ID i = 0; i < top; i++) cnt[i] = 0;
  for (ID e = 0; e < top - 1; e++) {
    const ID s = edges[e].src;
    const ID d = edges[e].dst;
    if ((x[d] != x[s]) || (y[d] != y[s])) {
      adj[s][cnt[s]] = e;
      cnt[s]++;
      adj[d][cnt[d]] = e;
      cnt[d]++;
    }
  }

  // find best distance for each triangle
  for (ID e = 0; e < top - 1; e++) dist[e] = -1;
  for (ID s = 0; s < top; s++) {
    if (cnt[s] >= 2) {
      const ctype x0 = x[s];
      const ctype y0 = y[s];
      for (char j = 0; j < cnt[s] - 1; j++) {
        const ID e1 = adj[s][j];
        const ID d1 = (s != edges[e1].src) ? edges[e1].src : edges[e1].dst;
        const ctype x1 = x[d1];
        const ctype y1 = y[d1];
        for (char k = j + 1; k < cnt[s]; k++) {
          const ID e2 = adj[s][k];
          const ID d2 = (s != edges[e2].src) ? edges[e2].src : edges[e2].dst;
          const ctype stx = std::max(std::min(x0, x1), std::min(std::max(x0, x1), x[d2]));
          const ctype sty = std::max(std::min(y0, y1), std::min(std::max(y0, y1), y[d2]));
          const ctype rd = std::abs(stx - x0) + std::abs(sty - y0);
          if (rd > 0) {
            const ctype rd1 = rd * (MaxPins * 2) + e1;  // tie breaker
            const ctype rd2 = rd * (MaxPins * 2) + e2;  // tie breaker
            if (dist[e1] < rd2) dist[e1] = rd2;
            if (dist[e2] < rd1) dist[e2] = rd1;
          }
        }
      }
    }
  }

  // process "triangles" to find best candidate Steiner points
  bool updated = false;
  for (ID e1 = 0; e1 < top - 2; e1++) {
    const ctype d1 = dist[e1];
    if (d1 > 0) {
      const ID e2 = d1 % (MaxPins * 2);
      if (e2 > e1) {
        const ctype d2 = dist[e2];
        if (e1 == d2 % (MaxPins * 2)) {
          const ctype x0 = x[edges[e1].src];
          const ctype y0 = y[edges[e1].src];
          const ctype x1 = x[edges[e1].dst];
          const ctype y1 = y[edges[e1].dst];
          ctype x2 = x[edges[e2].src];
          ctype y2 = y[edges[e2].src];
          if (((x2 == x0) && (y2 == y0)) || ((x2 == x1) && (y2 == y1))) {
            x2 = x[edges[e2].dst];
            y2 = y[edges[e2].dst];
          }
          updated = true;
          x[num] = std::max(std::min(x0, x1), std::min(std::max(x0, x1), x2));
          y[num] = std::max(std::min(y0, y1), std::min(std::max(y0, y1), y2));
          num++;
        }
      }
    }
  }

  return updated;
}


static inline void processNet(const int pins, const ctype* const __restrict__ xin, const ctype* const __restrict__ yin, ctype* const __restrict__ xout, ctype* const __restrict__ yout, edge* const __restrict__ edges)
{
  // initialize arrays and copy input coords to output
  for (ID j = 0; j < pins; j++) xout[j] = xin[j];
  for (ID j = 0; j < pins; j++) yout[j] = yin[j];

  // process nets
  if (pins == 2) {
    edges[0] = edge{0, 1};
  } else if (pins == 3) {
    ctype x0 = xout[0];
    ctype x1 = xout[1];
    ctype y0 = yout[0];
    ctype y1 = yout[1];
    xout[3] = std::max(std::min(x0, x1), std::min(std::max(x0, x1), xout[2]));
    yout[3] = std::max(std::min(y0, y1), std::min(std::max(y0, y1), yout[2]));
    edges[0] = edge{0, 3};
    edges[1] = edge{1, 3};
    edges[2] = edge{2, 3};
  } else {
    // iterate until all Steiner points added
    ID cnt = pins;
    int iter = 0;
    do {
      iter++;
      buildMST(cnt, xout, yout, edges);
    } while (insertSteinerPoints(cnt, xout, yout, edges));
    printf("%d iterations\n", iter);
  }
}


static void computeRSMT(const ctype* const __restrict__ xin, const ctype* const __restrict__ yin, ctype* const __restrict__ xout, ctype* const __restrict__ yout, edge* const __restrict__ edges, const int pins)
{
  // start time
  timeval start, end;
  gettimeofday(&start, NULL);

  // compute Steiner points and edges
  const int size = 2 * pins;
  for (int j = 0; j < size; j++) edges[j] = edge{0, 0};
  processNet(pins, xin, yin, xout, yout, edges);

  // end time
  gettimeofday(&end, NULL);
  const double runtime = end.tv_sec - start.tv_sec + (end.tv_usec - start.tv_usec) / 1000000.0;
  printf("compute time: %.6f s\n", runtime);
}


static ctype treeLength(const ID num, const ctype* const __restrict__ x, const ctype* const __restrict__ y, const edge* const __restrict__ edges)
{
  // compute wire length of Steiner tree
  ctype len = 0;
  for (ID i = 0; i < num - 1; i++) {
    const ctype x1 = x[edges[i].src];
    const ctype y1 = y[edges[i].src];
    const ctype x2 = x[edges[i].dst];
    const ctype y2 = y[edges[i].dst];
    len += std::abs(x1 - x2) + std::abs(y1 - y2);
  }
  return len;
}


// source of hash function: https://stackoverflow.com/questions/664014/what-integer-hash-function-are-good-that-accepts-an-integer-hash-key
static inline unsigned int hash(unsigned int val)
{
  val = ((val >> 16) ^ val) * 0x45d9f3b;
  val = ((val >> 16) ^ val) * 0x45d9f3b;
  return (val >> 16) ^ val;
}


int main(int argc, char* argv[])
{
  printf("\nSFP single net (%s)\n", __FILE__);
  printf("Copyright 2019-2022 Texas State University\n\n");

  // check command line
  if (argc != 2) {printf("USAGE: %s number_of_pins\n", argv[0]); exit(-1);}
  const int pins = atoi(argv[1]);
  if (pins > MaxPins) {printf("ERROR: too many pins\n"); exit(-1);}

  // random pin coordinates (change to actual coordinates)
  ctype* const xin = new ctype [pins];
  ctype* const yin = new ctype [pins];
  const int width = 1000;
  for (ID j = 0; j < pins; j++) xin[j] = hash(j + 1) % width;
  for (ID j = 0; j < pins; j++) yin[j] = hash(-j - 1) % width;

  // allocate result storage
  const int size = 2 * pins;
  ctype* const xout = new ctype [size];
  ctype* const yout = new ctype [size];
  edge* const edges = new edge [size];

  // compute Steiner points and edges
  computeRSMT(xin, yin, xout, yout, edges, pins);

  // print result
  const ctype len = treeLength(size, xout, yout, edges);  // body of treeLength function illustrates how to read solution
  printf("resulting tree length: %d\n", len);

  // clean up
  delete [] edges;  delete [] yout;  delete [] xout;
  delete [] yin;  delete [] xin;

  return 0;
}