/*
ECL MIS code:
This code finds the maximal independent set of an input. It uses the complex
encoding scheme discussed in the presentation. In this version of the code we
employ a modified encoding scheme to lower the number of bit-flips in the final
data structure

Copyright 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 and Alex Fallin

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

#include <stdlib.h>
#include <stdio.h>
#include <cuda.h>
#include <sys/time.h>
#include "ECLgraph.h"
#include "gpu_energy_monitor.h"

static const int Device = 0;
static const int ThreadsPerBlock = 256;

typedef unsigned char stattype;
static const stattype bfr_in = 0x80;
static const stattype out = 0;
static const stattype in = 0xfe;

static void CheckCuda() {
  cudaError_t e;
  cudaDeviceSynchronize();
  if (cudaSuccess != (e = cudaGetLastError())) {
    fprintf(stderr, "CUDA error %d: %s\n", e, cudaGetErrorString(e));
    exit(-1);
  }
}

/* main computation kernel */
static __global__ void findmins(const int nodes, const int* const __restrict__ nidx, const int* const __restrict__ nlist, volatile stattype* const __restrict__ nstat) {
  const int from = threadIdx.x + blockIdx.x * ThreadsPerBlock;
  const int incr = gridDim.x * ThreadsPerBlock;

  int missing;
  do {
    missing = 0;
    for (int v = from; v < nodes; v += incr) {
      const stattype nv = nstat[v];
      if (nv & 1) {
        int i = nidx[v];
        stattype temp;
        while ((i < nidx[v + 1]) && ((nv > (temp = nstat[nlist[i]])) || ((nv == temp) && (v > nlist[i])))) {
          i++;
        }
        if (i < nidx[v + 1]) {
          missing = 1;
        } else {
          for (int i = nidx[v]; i < nidx[v + 1]; i++) {
            nstat[nlist[i]] = out;
          }
          nstat[v] = in;
        }
      }
    }
  } while (missing != 0);
}

static __global__ void bfr_findmins(const int nodes, const int* const __restrict__ nidx, const int* const __restrict__ nlist, volatile stattype* const __restrict__ nstat) {
  const int from = threadIdx.x + blockIdx.x * ThreadsPerBlock;
  const int incr = gridDim.x * ThreadsPerBlock;

  int missing;
  do {
    missing = 0;
    for (int v = from; v < nodes; v += incr) {
      const stattype nv = nstat[v];
      if (nv & 0x7F) {
        int i = nidx[v];
        stattype temp;
        while ((i < nidx[v + 1]) && ((nv > (temp = nstat[nlist[i]])) || ((nv == temp) && (v > nlist[i])))) {
          i++;
        }
        if (i < nidx[v + 1]) {
          missing = 1;
        } else {
          for (int i = nidx[v]; i < nidx[v + 1]; i++) {
            nstat[nlist[i]] = out;
          }
          nstat[v] = bfr_in;
        }
      }
    }
  } while (missing != 0);
}

/* hash function to generate random values */

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

/* prioritized-selection initialization kernel */

static __global__ void init(const int nodes, const int edges, const int* const __restrict__ nidx, stattype* const __restrict__ nstat) {
  const int from = threadIdx.x + blockIdx.x * ThreadsPerBlock;
  const int incr = gridDim.x * ThreadsPerBlock;

  const float avg = (float) edges / nodes;
  const float scaledavg = ((in / 2) - 1) * avg;

  for (int i = from; i < nodes; i += incr) {
    stattype val = in;
    const int degree = nidx[i + 1] - nidx[i];
    if (degree > 0) {
      float x = degree - (hash(i) * 0.00000000023283064365386962890625f);
      int res = __float2int_rn(scaledavg / (avg + x));
      val = (res + res) | 1;
    }
    nstat[i] = val;
  }
}

static __global__ void bfr_init(const int nodes, const int edges, const int* const __restrict__ nidx, stattype* const __restrict__ nstat) {
  const int from = threadIdx.x + blockIdx.x * ThreadsPerBlock;
  const int incr = gridDim.x * ThreadsPerBlock;

  const float avg = (float) edges / nodes;
  const float scaledavg = ((in / 2) - 1) * avg;

  for (int i = from; i < nodes; i += incr) {
    stattype val = bfr_in;
    const int degree = nidx[i + 1] - nidx[i];
    if (degree > 0) {
      float x = degree - (hash(i) * 0.00000000023283064365386962890625f);
      int res = __float2int_rn(scaledavg / (avg + x));
      val = res + 1;
    }
    nstat[i] = val;
  }
}

struct GPUTimer {
  cudaEvent_t beg, end;
  GPUTimer() {
    cudaEventCreate(&beg);
    cudaEventCreate(&end);
  }
  ~GPUTimer() {
    cudaEventDestroy(beg);
    cudaEventDestroy(end);
  }
  void start() { cudaEventRecord(beg, 0); }
  float stop() {
    cudaEventRecord(end, 0);
    cudaEventSynchronize(end);
    float ms;
    cudaEventElapsedTime(&ms, beg, end);
    return 0.001f * ms;
  }
};

static void computeMIS(const int nodes, const int edges, const int* const __restrict__ nidx, const int* const __restrict__ nlist, stattype* const __restrict__ nstat) {
  cudaSetDevice(Device);
  cudaDeviceProp deviceProp;
  cudaGetDeviceProperties(&deviceProp, Device);
  if ((deviceProp.major == 9999) && (deviceProp.minor == 9999)) {
    fprintf(stderr, "ERROR: there is no CUDA capable device\n\n");
    exit(-1);
  }
  const int SMs = deviceProp.multiProcessorCount;
  const int mTpSM = deviceProp.maxThreadsPerMultiProcessor;
  printf("gpu: %s with %d SMs and %d mTpSM (%.1f MHz core and %.1f MHz mem)\n", deviceProp.name, SMs, mTpSM, deviceProp.clockRate * 0.001, deviceProp.memoryClockRate * 0.001);

  int* nidx_d;
  int* nlist_d;
  stattype* nstat_d;

  if (cudaSuccess != cudaMalloc((void**) &nidx_d, (nodes + 1) * sizeof(int))) {
    fprintf(stderr, "ERROR: could not allocate nidx_d\n\n");
    exit(-1);
  }
  if (cudaSuccess != cudaMalloc((void**) &nlist_d, edges * sizeof(int))) {
    fprintf(stderr, "ERROR: could not allocate nlist_d\n\n");
    exit(-1);
  }
  if (cudaSuccess != cudaMalloc((void**) &nstat_d, nodes * sizeof(stattype))) {
    fprintf(stderr, "ERROR: could not allocate nstat_d\n\n");
    exit(-1);
  }

  if (cudaSuccess != cudaMemcpy(nidx_d, nidx, (nodes + 1) * sizeof(int), cudaMemcpyHostToDevice)) {
    fprintf(stderr, "ERROR: copying to device failed\n\n");
    exit(-1);
  }
  if (cudaSuccess != cudaMemcpy(nlist_d, nlist, edges * sizeof(int), cudaMemcpyHostToDevice)) {
    fprintf(stderr, "ERROR: copying to device failed\n\n");
    exit(-1);
  }

  cudaFuncSetCacheConfig(init, cudaFuncCachePreferL1);
  cudaFuncSetCacheConfig(findmins, cudaFuncCachePreferL1);

  const int blocks = SMs * mTpSM / ThreadsPerBlock;

  bool nstat_reg_f = false;
  stattype* nstat_reg = (stattype*) malloc(nodes * sizeof(nstat[0]));

  // Running all kernels to get it all CUDA initialized
  init<<<blocks, ThreadsPerBlock>>>(nodes, edges, nidx_d, nstat_d);
  findmins<<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nlist_d, nstat_d);
  cudaDeviceSynchronize();
  if (cudaSuccess != cudaMemcpy(nstat, nstat_d, nodes * sizeof(stattype), cudaMemcpyDeviceToHost)) {
    fprintf(stderr, "ERROR: copying from device failed\n\n");
    exit(-1);
  }
  bfr_init<<<blocks, ThreadsPerBlock>>>(nodes, edges, nidx_d, nstat_d);
  bfr_findmins<<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nlist_d, nstat_d);
  cudaDeviceSynchronize();

  // Finding number of reps
  timeval beg, end;
  int reps = 0;
  gettimeofday(&beg, NULL);
  do {
    init<<<blocks, ThreadsPerBlock>>>(nodes, edges, nidx_d, nstat_d);
    findmins<<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nlist_d, nstat_d);
    cudaDeviceSynchronize();

    gettimeofday(&end, NULL);
    reps++;
  } while ((end.tv_sec - beg.tv_sec + (end.tv_usec - beg.tv_usec) / 1000000.0) < 10.0);
  CheckCuda();

  // Run experiments interleaved 10 times
  const int tot_runs = 10;

  // GPU Energy Setup
  GPUMonitor gpu;
  gpu.initMonitor(Device);

  // Trackers
  unsigned long long regular_consumption[tot_runs];
  double regular_runtime[tot_runs];
  unsigned long long bfr_consumption[tot_runs];
  double bfr_runtime[tot_runs];

  // Perform experiments
  for (int run = 0; run < tot_runs; run++) {
    // Regular
    gpu.startEnergy();
    gettimeofday(&beg, NULL);
    for (int rep = 0; rep < reps; rep++) {
      init<<<blocks, ThreadsPerBlock>>>(nodes, edges, nidx_d, nstat_d);
      findmins<<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nlist_d, nstat_d);
      cudaDeviceSynchronize();
    }
    gettimeofday(&end, NULL);
    regular_consumption[run] = gpu.getEnergyConsumption();
    regular_runtime[run] = end.tv_sec - beg.tv_sec + (end.tv_usec - beg.tv_usec) / 1000000.0;
    if (!nstat_reg_f){
      nstat_reg_f = true;
      if (cudaSuccess != cudaMemcpy(nstat_reg, nstat_d, nodes * sizeof(stattype), cudaMemcpyDeviceToHost)) {fprintf(stderr, "ERROR: copying from device failed\n\n");  exit(-1);}
    }

    // BFR
    gpu.startEnergy();
    gettimeofday(&beg, NULL);
    for (int rep = 0; rep < reps; rep++) {
      bfr_init<<<blocks, ThreadsPerBlock>>>(nodes, edges, nidx_d, nstat_d);
      bfr_findmins<<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nlist_d, nstat_d);
      cudaDeviceSynchronize();
    }
    gettimeofday(&end, NULL);
    bfr_consumption[run] = gpu.getEnergyConsumption();
    bfr_runtime[run] = end.tv_sec - beg.tv_sec + (end.tv_usec - beg.tv_usec) / 1000000.0;
  }

  // Print results
  // reps, reg_con, reg_time, bfr_con, bfr_time
  for (int run = 0; run < tot_runs; run++) {
    printf("%d, %lld, %.4f, %lld, %.4f\n", reps, regular_consumption[run], regular_runtime[run], bfr_consumption[run], bfr_runtime[run]);
  }

  if (cudaSuccess != cudaMemcpy(nstat, nstat_d, nodes * sizeof(stattype), cudaMemcpyDeviceToHost)) {fprintf(stderr, "ERROR: copying from device failed\n\n");  exit(-1);}

  for (int i = 0; i < nodes; i++) {
    if ((nstat_reg[i] == in) != (nstat[i] == bfr_in)) {
      printf("Missmatch at %d\n", i);
    }
  }

  delete [] nstat_reg;
  cudaFree(nstat_d);
  cudaFree(nlist_d);
  cudaFree(nidx_d);
}

int main(int argc, char* argv[]) {
  printf("ECL-MIS v1.3 (%s)\n", __FILE__);
  printf("Copyright 2017-2022 Texas State University\n");

  if (argc != 2) {
    fprintf(stderr, "USAGE: %s input_file_name\n\n", argv[0]);
    exit(-1);
  }

  ECLgraph g = readECLgraph(argv[1]);
  printf("%s\n", argv[1]);
  printf("configuration: %d nodes and %d edges (%s)\n", g.nodes, g.edges, argv[1]);
  printf("average degree: %.2f edges per node\n", 1.0 * g.edges / g.nodes);

  stattype* nstatus = (stattype*) malloc(g.nodes * sizeof(nstatus[0]));
  if (nstatus == NULL) {
    fprintf(stderr, "ERROR: could not allocate nstatus\n\n");
    exit(-1);
  }

  computeMIS(g.nodes, g.edges, g.nindex, g.nlist, nstatus);

  int count = 0;
  for (int v = 0; v < g.nodes; v++) {
    if (nstatus[v] == bfr_in) {
      count++;
    }
  }
  printf("elements in set: %d (%.1f%%)\n", count, 100.0 * count / g.nodes);

  freeECLgraph(g);
  free(nstatus);
  return 0;
}