/*
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 = 256;  // 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 inline void buildMST(const ID num, const ctype* const __restrict__ x, const ctype* const __restrict__ y, edge* const __restrict__ edges)
{
  ctype dist[MaxPins * 2];
  ID source[MaxPins * 2];
  ID destin[MaxPins * 2];

  // 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)
{
  ctype dist[MaxPins * 2];
  ID adj[MaxPins * 2][8];
  char cnt[MaxPins * 2];
  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 i, const int* const __restrict__ idxin, const ctype* const __restrict__ xin, const ctype* const __restrict__ yin, int* const __restrict__ idxout, ctype* const __restrict__ xout, ctype* const __restrict__ yout, edge* const __restrict__ edges)
{
  // initialize arrays and copy input coords to output
  const int pin = idxin[i];
  const ID num = idxin[i + 1] - pin;
  const int pout = 2 * pin;
  idxout[i] = pout;
  for (ID j = 0; j < num; j++) xout[pout + j] = xin[pin + j];
  for (ID j = 0; j < num; j++) yout[pout + j] = yin[pin + j];

  // process nets
  if (num == 2) {
    edges[pout] = edge{0, 1};
  } else if (num == 3) {
    ctype x0 = xout[pout];
    ctype x1 = xout[pout + 1];
    ctype y0 = yout[pout];
    ctype y1 = yout[pout + 1];
    xout[pout + 3] = std::max(std::min(x0, x1), std::min(std::max(x0, x1), xout[pout + 2]));
    yout[pout + 3] = std::max(std::min(y0, y1), std::min(std::max(y0, y1), yout[pout + 2]));
    edges[pout] = edge{0, 3};
    edges[pout + 1] = edge{1, 3};
    edges[pout + 2] = edge{2, 3};
  } else {
    // iterate until all Steiner points added
    ID cnt = num;
    do {
      buildMST(cnt, &xout[pout], &yout[pout], &edges[pout]);
    } while (insertSteinerPoints(cnt, &xout[pout], &yout[pout], &edges[pout]));
  }
}


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

  // compute Steiner points and edges
  const int size = 2 * idxin[numnets];
  idxout[numnets] = size;

  #pragma omp parallel
  {
    #pragma omp for
    for (int j = 0; j < size; j++) edges[j] = edge{0, 0};

    // process nets
    #pragma omp for schedule(dynamic, 256)  // to avoid false sharing and WL overhead
    for (int i = 0; i < numnets; i++) {
      processNet(i, idxin, xin, yin, idxout, 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);
  printf("throughput: %.f nets/sec\n", numnets / 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;
}


// struct to store grid information
struct grid {
  int grid[3];
  int* vc = NULL;
  int* hc = NULL;
  int* min_wid = NULL;
  int* min_space = NULL;
  int* via_space = NULL;
  int llx, lly, tile_wid, tile_height;
};


// struct to store net and pin information
struct net_list
{
  int num_net;
  std::vector<std::tuple<int, int> >* num_net_arr = NULL;
  int* net_id = NULL;
  int* net_num_pins = NULL;
  int* net_min_wid = NULL;
};


// function to read in input file
static void read_file(const char* file, grid& g, net_list& n)
{
  std::string line;
  std::string text1, text2;
  int line_count = 0;  // for error messages

  std::fstream myfile(file);
  if (!myfile.is_open()) {std::cout << "ERROR: Cannot open input file!\n"; exit(-1);}

  // read grid x and y co-ordinates and number of layers
  getline(myfile, line);
  line_count++;
  std::stringstream data1(line);
  if (!(data1 >> text1 >> g.grid[0] >> g.grid[1] >> g.grid[2])) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if (text1 != "grid") {std::cout << "ERROR: Invalid format of grid!\n"; exit(-1);}
  for (int i = 0; i < 3; i++) {
    if (g.grid[i] < 1) {std::cout << "ERROR: Grid data should be a reasonable number!\n"; exit(-1);}
  }

  // read vertical capacity of each layer
  getline(myfile, line);
  line_count++;
  std::stringstream data2(line);
  g.vc = new int [g.grid[2] + 1];
  if (!(data2 >> text1 >> text2)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if ((text1 != "vertical") || (text2 != "capacity")) {std::cout << "ERROR: Invalid format of g.vertical capacity!\n"; exit(-1);}
  for (int i = 1; i <= g.grid[2]; i++) {
      if (data2 >> g.vc[i]) {
      if (g.vc[i] < 0) {std::cout << "ERROR: vertical capacity should be a reasonable number!\n"; exit(-1);}
    }
  }

  // read horizontal capacity of each layer
  getline(myfile, line);
  line_count++;
  std::stringstream data3(line);
  g.hc = new int [g.grid[2] + 1];
  if (!(data3 >> text1 >> text2)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n";}
  if ((text1 != "horizontal") || (text2 != "capacity")) {std::cout << "ERROR: Invalid format of g.horizontal capacity!\n"; exit(-1);}
  for (int i = 1; i <= g.grid[2]; i++) {
    if (data3 >> g.hc[i]) {
      if (g.hc[i] < 0) {std::cout << "ERROR: horizontal capacity should be a reasonable number!\n"; exit(-1);}
    }
  }

  // read minimum width of each layer
  getline(myfile, line);
  line_count++;
  std::stringstream data4(line);
  g.min_wid = new int [g.grid[2] + 1];
  if (!(data4 >> text1 >> text2)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if ((text1 != "minimum") || (text2 != "width")) {std::cout << "ERROR: Invalid format of minimum width!\n"; exit(-1);}
  for (int i = 1; i <= g.grid[2] + 1; i++) {
    if (data4 >> g.min_wid[i]) {
      if (g.min_wid[i] < 1) {std::cout << "ERROR: Minimum width should be a reasonable number!\n"; exit(-1);}
    }
  }

  // read minimum spacing of each layer
  getline(myfile, line);
  line_count++;
  std::stringstream data5(line);
  g.min_space = new int [g.grid[2] + 1];
  if (!(data5 >> text1 >> text2)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if ((text1 != "minimum") || (text2 != "spacing")) {std::cout << "ERROR: Invalid format of minimum spacing!\n"; exit(-1);}
  for (int i = 1; i <= g.grid[2]; i++) {
    if (data5 >> g.min_space[i]) {
      if (g.min_space[i] < 0) {std::cout << "ERROR: Minimum spacing should be a reasonable number!\n"; exit(-1);}
    }
  }

  // read via spacing of each layer
  getline(myfile, line);
  line_count++;
  std::stringstream data6(line);
  g.via_space = new int [g.grid[2] + 1];
  if (!(data6 >> text1 >> text2)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if ((text1 != "via") || (text2 != "spacing")) {std::cout << "ERROR: Invalid format of via spacing!\n"; exit(-1);}
  for (int i = 1; i <= g.grid[2]; i++) {
    if (data6 >> g.via_space[i]) {
      if (g.via_space[i] < 0) {std::cout << "ERROR: Via spacing should be a reasonable number!\n"; exit(-1);}
    }
  }

  // read lower left x and y co-ordinates for the global routing region, tile width and tile height per layer
  getline(myfile, line);
  line_count++;
  std::stringstream data7(line);
  if (!(data7 >> g.llx >> g.lly >> g.tile_wid >> g.tile_height)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if (g.tile_wid < 1 || g.tile_height < 1) {std::cout << "ERROR: Tile width and tile height should be a reasonable number!\n"; exit(-1);}

  // read total number of nets
  do {getline(myfile, line);} while (line == "");
  line_count++;
  std::stringstream data8(line);
  if (!(data8 >> text1 >> text2 >> n.num_net)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
  if ((text1 != "num") || (text2 != "net")) {std::cout << "ERROR: Invalid format of num net!\n"; exit(-1);}
  if (n.num_net < 1) {std::cout << "ERROR: Number of nets should be a reasonable number!\n";}

  // allocate memory
  n.net_id = new int [n.num_net];
  n.net_num_pins = new int [n.num_net];
  n.net_min_wid = new int [n.num_net];
  n.num_net_arr = new std::vector<std::tuple<int, int> > [n.num_net];

  // read net name, net id, number of pins and minimum width of each net
  for (int i = 0; i < n.num_net; i++) {
    getline(myfile, line);
    line_count++;
    std::stringstream data9(line);
    if (!(data9 >> text1 >> n.net_id[i] >> n.net_num_pins[i] >> n.net_min_wid[i])) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
    if ((n.net_id[i] < 0) || (n.net_num_pins[i] < 2) || (n.net_min_wid[i] < 1)) {std::cout << "ERROR: Net ID, number of pins and min width should be a reasonable number!\n"; exit(-1);}

    // read x, y and layer information of each net
    for (int j = 0; j < n.net_num_pins[i]; j++) {
      getline(myfile, line);
      line_count++;
      std::stringstream data10(line);
      int x, y, layer;
      if (!(data10 >> x >> y >> layer)) {std::cout << "ERROR: Line" << line_count << " couldn't be parsed!\n"; exit(-1);}
      if (!(x >= g.llx && x < (g.grid[0] * g.tile_wid))) x = (g.grid[0] * g.tile_wid) - 1;
      if (!(y >= g.lly && y < (g.grid[1] * g.tile_height))) y = (g.grid[1] * g.tile_height) - 1;
      if (!(layer > 0 && layer <= g.grid[2])) {std::cout << "ERROR: layer should be within grid!\n"; exit(-1);}
      if ((x < 0) || (y < 0) || (x >= INT_MAX / (MaxPins * 4)) || (y >= INT_MAX / (MaxPins * 4))) {std::cout << "ERROR: x or y out of bounds\n"; exit(-1);}
      n.num_net_arr[i].push_back(std::make_tuple(x, y));
    }
  }

  myfile.close();
}


// free info from input file
static void free_memory(const grid& g, const net_list& n)
{
  delete [] g.vc;
  delete [] g.hc;
  delete [] g.min_wid;
  delete [] g.min_space;
  delete [] g.via_space;
  delete [] n.net_id;
  delete [] n.net_num_pins;
  delete [] n.net_min_wid;
  delete [] n.num_net_arr;
}


int main(int argc, char* argv[])
{
  #ifdef _OPENMP
    printf("\nSFP OpenMP (%s)\n", __FILE__);
  #else
    printf("\nSFP serial (%s)\n", __FILE__);
  #endif
  printf("Copyright 2019-2022 Texas State University\n\n");

  // check command line
  if (argc != 3) {printf("USAGE: %s file_name print_results\n", argv[0]); exit(-1);}
  const bool print = (atoi(argv[2]) != 0);

  // read input file
  printf("reading: %s\n", argv[1]);
  grid g;
  net_list n;
  read_file(argv[1], g, n);
  const int numnets = n.num_net;

  // start converting data structure
  int* const idxin = new int [numnets + 1];
  idxin[0] = 0;
  ID hipin = 0;
  int pos = 0;
  for (int i = 0; i < numnets; i++) {
    const ID num = std::min(n.net_num_pins[i], MaxPins);
    hipin = std::max(hipin, num);
    pos += num;
    idxin[i + 1] = pos;
  }

  // print histogram
  int trunc = 0;
  int* const hist = new int [hipin + 1];
  for (int i = 0; i < hipin + 1; i++) hist[i] = 0;
  for (int i = 0; i < numnets; i++) {
    const ID num = std::min(n.net_num_pins[i], MaxPins);
    if (num < n.net_num_pins[i]) trunc++;
    hist[num]++;
  }
  int sum = 0;
  for (int i = 0; i < hipin + 1; i++) {
    sum += hist[i];
    if (print) printf("hist %4d: %6.2f%%  %6.2f%%  %d\n", i, 100.0 * hist[i] / numnets, 100.0 * sum / numnets, hist[i]);
  }
  delete [] hist;

  // print info
  if (print) printf("\n");
  printf("number of nets: %d\n", numnets);
  printf("max pins per net: %d\n", hipin);
  printf("truncated nets: %d\n", trunc);
  if (hipin > MaxPins) {printf("ERROR: hi_pin_count must be no more than %d\n", MaxPins); exit(-1);}

  // copy pin coordinates
  ctype* const xin = new ctype [idxin[numnets]];
  ctype* const yin = new ctype [idxin[numnets]];
  pos = 0;
  for (int i = 0; i < numnets; i++) {
    const ID num = idxin[i + 1] - idxin[i];
    for (ID j = 0; j < num; j++) xin[pos + j] = std::get<0>(n.num_net_arr[i][j]);
    for (ID j = 0; j < num; j++) yin[pos + j] = std::get<1>(n.num_net_arr[i][j]);
    pos += num;
  }

  // allocate result storage
  const int size = 2 * idxin[numnets];
  int* const idxout = new int [numnets + 1];
  ctype* const xout = new ctype [size];
  ctype* const yout = new ctype [size];
  edge* const edges = new edge [size];

  // compute Steiner points and edges
  computeRSMT(idxin, xin, yin, idxout, xout, yout, edges, numnets);

  // print result
  long total = 0;
  if (print) printf("\nresulting tree lengths:\n");
  for (int i = 0; i < numnets; i++) {
    const ctype len = treeLength(idxout[i + 1] - idxout[i], &xout[idxout[i]], &yout[idxout[i]], &edges[idxout[i]]);  // body of treeLength function illustrates how to read solution
    total += len;
    if (print) {
      const ID num = std::min(n.net_num_pins[i], MaxPins);
      printf("%d: %d\n", num, len);
    }
  }
  printf("\ntotal wirelength: %ld\n", total);

  // clean up
  free_memory(g, n);
  delete [] edges;  delete [] yout;  delete [] xout;  delete [] idxout;
  delete [] yin;  delete [] xin;  delete [] idxin;

  return 0;
}