#include <map>
#include <deque>
#include <vector>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <sys/time.h>

#include <cmath>
#include <cstring>

#ifdef __APPLE_CC__
#include <OpenCL/opencl.h>
#else
#include <CL/opencl.h>
#endif

using namespace std;

namespace xppc{
#define OFLA         // omit the flasher DOM
#define ROMB         // use rhomb cells aligned with the array

#define ASENS        // enable angular sensitivity
#define TILT         // enable tilted ice layers
//#define RAND         // disable for deterministic results

#define MKOW         // use Marek Kowalski's photon yield parametrization
#define ANGW         // smear cherenkov cone due to shower development
#define LONG         // simulate longitudinal cascade development

#ifdef ASENS
#define ANUM 11      // number of coefficients in the angular sensitivity curve
#endif

#ifdef TILT
#define LMAX 6       // number of dust loggers
#define LYRS 170     // number of depth points
#endif

#ifdef ROMB
#define DIR1 9.3
#define DIR2 129.3

#define CX   21
#define CY   19
#define NSTR 94
#else
#define CX   13
#define CY   13
#define NSTR 94
#endif

#define OVER 10      // size of photon bunches along the muon track

#define HQUO   16    // save at most photons/HQUO hits
#define NPHO   1024  // maximum number of photons propagated by one thread

#define WNUM   32    // number of wavelength slices
#define MAXLYS 171   // maximum number of ice layers
#define MAXGEO 5200  // maximum number of OMs
#define MAXRND 16028 // maximum number of random number multipliers

#define XXX 1.e-5f
#define FPI 3.141592653589793f
#define OMR 0.16510f // DOM radius [m]

#include "pro.cxx"
#include "ini.cxx"

#define SHRT
#include "pro.cxx"
#undef SHRT

#ifdef XLIB
  void igeo(float enh){ m.geo(); d.eff*=enh; }
#endif

  unsigned int pmax, pmxo, pn;

  int nblk, nthr;
  void ini(){
    pmax=nblk*nthr*NPHO;
    d.hnum=pmax/HQUO;
    pmxo=pmax/OVER;
    pn=0;
  }

  float deviceTime=0;

  cl_device_id devID;
  cl_context ctx;
  cl_command_queue cq;
  cl_program program;
  cl_kernel clkernel;

  cl_mem ed, ez, eh, ep, eo; // pointers to structures on device

  void checkError(cl_int result){
    if(result!=CL_SUCCESS){
      cerr<<"OpenCL Error: "<<result<<endl;
      exit(2);
    }
  }

  void ini(int type){
    ini();

    {
      unsigned int size=d.rsize, need=nblk*nthr;
      if(size<need) cerr<<"Error: not enough multipliers: only have "<<size<<" (need "<<need<<")!"<<endl;
      else d.rsize=need;
    }

    int tot=0, cnt=0;
    cl_int err;

    {
      unsigned int size=sizeof(datz); tot+=size;
      ez = clCreateBuffer(ctx, CL_MEM_READ_WRITE|CL_MEM_COPY_HOST_PTR, size, &z, &err); checkError(err);
    }

    {
      q.hits = new hit[d.hnum];
      unsigned int size=d.hnum*sizeof(hit); tot+=size;
      eh = clCreateBuffer(ctx, CL_MEM_READ_WRITE, size, NULL, &err); checkError(err);
    }

    {
      unsigned int size=sizeof(photon);
      if(d.type==0){
	q.pz = new photon[pmxo];
	size*=pmxo; tot+=size;
      }
      ep = clCreateBuffer(ctx, CL_MEM_READ_WRITE, size, NULL, &err); checkError(err);
    }

    {
      unsigned int size=d.gsize*sizeof(DOM); cnt+=size;
      eo = clCreateBuffer(ctx, CL_MEM_READ_ONLY|CL_MEM_COPY_HOST_PTR, size, q.oms, &err); checkError(err);
    }

    {
      unsigned int size=sizeof(dats); tot+=size;
      ed = clCreateBuffer(ctx, CL_MEM_READ_WRITE|CL_MEM_COPY_HOST_PTR, size, &d, &err); checkError(err);
    }

    cerr<<"Total GPU memory usage: "<<tot<<"  const: "<<cnt<<endl;
  }

  void fin(){
    checkError(clReleaseMemObject(ez));
    checkError(clReleaseMemObject(eh));
    checkError(clReleaseMemObject(ep));
    checkError(clReleaseMemObject(eo));
    checkError(clReleaseMemObject(ed));

    if(d.type==0) delete q.pz;
    delete q.hits;
  }

  void listDevices(int device = -1){
    cl_uint num;
    checkError(clGetPlatformIDs(0, NULL, &num));
    fprintf(stderr, "\nFound %d platforms:\n", num);

    cl_platform_id ids[num];
    checkError(clGetPlatformIDs(num, ids, NULL));

    for(int i=0, n=0; i<num; ++i){
      size_t siz;
      cl_platform_id & platform=ids[i];
      fprintf(stderr, "platform %d: ", i);

      {
	checkError(clGetPlatformInfo(platform, CL_PLATFORM_NAME, 0, NULL, &siz));
	char buf[siz];
	checkError(clGetPlatformInfo(platform, CL_PLATFORM_NAME, siz, &buf, NULL));
	fprintf(stderr, "%s ", buf);
      }

      {
	checkError(clGetPlatformInfo(platform, CL_PLATFORM_VERSION, 0, NULL, &siz));
	char buf[siz];
	checkError(clGetPlatformInfo(platform, CL_PLATFORM_VERSION, siz, &buf, NULL));
	fprintf(stderr, "%s ", buf);
      }

      cl_uint num;
      checkError(clGetDeviceIDs(platform, CL_DEVICE_TYPE_ALL, 0, NULL, &num));
      fprintf(stderr, "device count: %d\n", num);

      cl_device_id devices[num];
      checkError(clGetDeviceIDs(platform, CL_DEVICE_TYPE_ALL, num, devices, NULL));

      for(int j=0; j<num; j++, n++){
	cl_device_id & di = devices[i];
	fprintf(stderr, "  device %d:", j);

	{
	  clGetDeviceInfo(di, CL_DEVICE_MAX_COMPUTE_UNITS, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_MAX_COMPUTE_UNITS, siz, buf, NULL);
	  fprintf(stderr, " cu=%d", * (cl_uint *) buf);
	}

	{
	  clGetDeviceInfo(di, CL_DEVICE_MAX_WORK_GROUP_SIZE, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_MAX_WORK_GROUP_SIZE, siz, buf, NULL);
	  fprintf(stderr, " gr=%lu", * (size_t *) buf);
	}

	{
	  clGetDeviceInfo(di, CL_DEVICE_MAX_CLOCK_FREQUENCY, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_MAX_CLOCK_FREQUENCY, siz, buf, NULL);
	  fprintf(stderr, " %d MHz", * (cl_uint *) buf);
	}

	{
	  clGetDeviceInfo(di, CL_DEVICE_GLOBAL_MEM_SIZE, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_GLOBAL_MEM_SIZE, siz, buf, NULL);
	  fprintf(stderr, " %lu bytes", * (cl_ulong *) buf);
	}

	{
	  clGetDeviceInfo(di, CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE, siz, buf, NULL);
	  fprintf(stderr, " cm=%lu", * (cl_ulong *) buf);
	}

	{
	  clGetDeviceInfo(di, CL_DEVICE_LOCAL_MEM_SIZE, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_LOCAL_MEM_SIZE, siz, buf, NULL);
	  fprintf(stderr, " lm=%lu", * (cl_ulong *) buf);
	}

	{
	  clGetDeviceInfo(di, CL_DEVICE_ERROR_CORRECTION_SUPPORT, 0, NULL, &siz);
	  char buf[siz];
	  clGetDeviceInfo(di, CL_DEVICE_ERROR_CORRECTION_SUPPORT, siz, buf, NULL);
	  fprintf(stderr, " ecc=%d", * (cl_bool *) buf);
	}

	if(n==device){
	  devID = di;
	  fprintf(stderr, " *");
	}

	fprintf(stderr, "\n");
      }
    }
  }

  static unsigned int old=0;

#include "f2k.cxx"

  void kernel(unsigned int num){
    if(old>0){
      checkError(clEnqueueReadBuffer(cq, ed, CL_TRUE, 0, 2*sizeof(int), &d, 0, NULL, NULL));
      if(d.ab>0) cerr<<"Error: TOT was a nan or an inf "<<d.ab<<" times!"<<endl;

      cerr<<"photons: "<<old<<"  hits: "<<d.hidx<<endl;
      if(d.hidx>=d.hnum){ d.hidx=d.hnum; cerr<<"Error: data buffer overflow occurred!"<<endl; }

      if(d.hidx>0){
	unsigned int size=d.hidx*sizeof(hit);
	checkError(clEnqueueReadBuffer(cq, eh, CL_TRUE, 0, size, q.hits, 0, NULL, NULL));
      }
    }

    if(num>0){
      unsigned int zero=0;
      checkError(clSetKernelArg(clkernel, 0, sizeof(uint), &zero));
      checkError(clSetKernelArg(clkernel, 1, sizeof(cl_mem), &ed));
      checkError(clSetKernelArg(clkernel, 2, sizeof(cl_mem), &ez));
      checkError(clSetKernelArg(clkernel, 3, sizeof(cl_mem), &eh));
      checkError(clSetKernelArg(clkernel, 4, sizeof(cl_mem), &ep));
      checkError(clSetKernelArg(clkernel, 5, sizeof(cl_mem), &eo));

      size_t ntot=1, nloc=1;
      checkError(clEnqueueNDRangeKernel(cq, clkernel, 1, NULL, &ntot, &nloc, 0, NULL, NULL));
      checkError(clFinish(cq));

      checkError(clSetKernelArg(clkernel, 0, sizeof(uint), &num));
      checkError(clSetKernelArg(clkernel, 1, sizeof(cl_mem), &ed));
      checkError(clSetKernelArg(clkernel, 2, sizeof(cl_mem), &ez));
      checkError(clSetKernelArg(clkernel, 3, sizeof(cl_mem), &eh));
      checkError(clSetKernelArg(clkernel, 4, sizeof(cl_mem), &ep));
      checkError(clSetKernelArg(clkernel, 5, sizeof(cl_mem), &eo));

      ntot=nblk*nthr, nloc=nthr;
      checkError(clEnqueueNDRangeKernel(cq, clkernel, 1, NULL, &ntot, &nloc, 0, NULL, NULL));
      checkError(clFinish(cq));
    }

    if(old>0) print();
    old=num;
  }

  void flini(int str, int dom){
    int type=1;
    float r[3]={0, 0, 0};

    if(str<0){ type=2; str=-str; }
    if(str==0) switch(dom){
      case 1: type=3; r[0]=544.07; r[1]=55.89; r[2]=136.86; break;
      case 2: type=4; r[0]=11.87; r[1]=179.19; r[2]=-205.64; break;
    }
    else for(int n=0; n<d.gsize; n++) if(q.names[n].str==str && q.names[n].dom==dom){
#ifdef OFLA
      d.fla=n;
#endif
      for(int m=0; m<3; m++) r[m]=q.oms[n].r[m]; break;
    }

    d.type=type; for(int m=0; m<3; m++) d.r[m]=r[m];
    ini(type);
  }

  void flone(unsigned long long num){
    for(long long i=num; i>0; i-=pmax) kernel(min(i, (long long) pmax));
    kernel(0);
  }

  void flasher(int str, int dom, unsigned long long num, int itr){
    flini(str, dom);

    for(int j=0; j<max(1, itr); j++){
      flone(num);
      if(itr>0) printf("\n");
    }

    fin();
  }

  void output(){
    if(pn>0){
      unsigned int size=pn*sizeof(photon);
      checkError(clEnqueueWriteBuffer(cq, ep, CL_TRUE, 0, size, q.pz, 0, NULL, NULL));
    }

    kernel(pn*OVER); pn=0;

    flnd=flne;
  }

  void start(){
  }

  void stop(){
    fprintf(stderr, "\nDevice time: %2.1f [ms]\n", deviceTime);
    if(clkernel) checkError(clReleaseKernel(clkernel));
    if(program) checkError(clReleaseProgram(program));
    if(cq) checkError(clReleaseCommandQueue(cq));
    if(ctx) checkError(clReleaseContext(ctx));
  }

  void choose(int device){
    sv+=device;
    listDevices(device);

    {
      bool verbose = false;

      cl_int err;

      ctx = clCreateContext(0, 1, &devID, NULL, NULL, &err); checkError(err);
      cq = clCreateCommandQueue(ctx, devID, 0, &err); checkError(err);

      string source;
      source.append("#pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable\n");
      const string tmp = kernel_source.substr(1, kernel_source.length()-2);

      if(verbose){
	int k=0, m=0;
	while((m=tmp.find("; ", k, 2))!=-1){
	  source.append(tmp, k, m-k+2);
	  source.append("\n");
	  k=m+1;
	}
	source.append(tmp, k, tmp.length());
      }
      else source.append(tmp);

      const char *src = source.c_str();
      size_t length=source.length();

      if(verbose) fprintf(stderr, "KERNEL SOURCE CODE:\n%s\n", src);
      program = clCreateProgramWithSource(ctx, 1, &src, &length, &err); checkError(err);

      clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
      checkError(clUnloadCompiler());

      cl_build_status status;
      checkError(clGetProgramBuildInfo(program, devID, CL_PROGRAM_BUILD_STATUS, sizeof(cl_build_status), &status, NULL));

      size_t siz;
      checkError(clGetProgramBuildInfo(program, devID, CL_PROGRAM_BUILD_LOG, 0, NULL, &siz));
      char log[siz+1]; log[siz] = '\0';
      checkError(clGetProgramBuildInfo(program, devID, CL_PROGRAM_BUILD_LOG, siz, log, NULL));
      if(verbose || status!=CL_SUCCESS) fprintf(stderr, "BUILD LOG:\n%s\n", log);
      checkError(status);

      clkernel = clCreateKernel(program, "propagate", &err); checkError(err);
    }

    {
      cl_uint count;
      clGetDeviceInfo(devID, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(cl_int), &count, NULL);
      nblk=count;

      size_t size;
      clGetKernelWorkGroupInfo(clkernel, devID, CL_KERNEL_WORK_GROUP_SIZE, sizeof(size_t), &size, NULL);
      nthr=size;

      cl_ulong lmem;
      clGetKernelWorkGroupInfo(clkernel, devID, CL_KERNEL_LOCAL_MEM_SIZE, sizeof(cl_ulong), &lmem, NULL);

      cerr<<"Running on "<<nblk<<" MPs x "<<nthr<<" threads; Kernel uses: l="<<lmem<<endl;
    }
  }
}

#ifndef XLIB
using namespace xppc;

int main(int arg_c, char *arg_a[]){
  start();
  if(arg_c<=1){
    listDevices();
    fprintf(stderr, "\nUse: %s [device] (f2k muons)\n"
	    "     %s [str] [om] [num] [device] (flasher)\n", arg_a[0], arg_a[0]);
  }
  else if(arg_c<=2){
    int device=0;
    if(arg_c>1) device=atoi(arg_a[1]);
    choose(device);
    fprintf(stderr, "Processing f2k muons from stdin on device %d\n", device);
    f2k();
  }
  else{
    int str=0, dom=0, device=0, itr=0;
    unsigned long long num=1000000ULL;

    if(arg_c>1) str=atoi(arg_a[1]);
    if(arg_c>2) dom=atoi(arg_a[2]);
    if(arg_c>3){
      num=(unsigned long long) atof(arg_a[3]);
      char * sub = strchr(arg_a[3], '*');
      if(sub!=NULL) itr=(int) atof(++sub);
    }
    if(arg_c>4) device=atoi(arg_a[4]);
    choose(device);
    fprintf(stderr, "Running flasher simulation on device %d\n", device);
    flasher(str, dom, (unsigned long long)llroundf(num*(long double)d.eff), itr);
  }

  stop();
}
#endif