#if defined __GNUC__ && !defined __INTEL_COMPILER
#include <omp.h>
#endif


#if defined __GNUC__ || defined __INTEL_COMPILER
#include <random>
#endif
#include <ctime>
#include <cstdlib>
#endif
#include "GMS_grass_scatterers_AVX.h"
#if (SAMPLE_HW_PMC) == 1
    #include "libpfc.h"
    #include <string.h>
    #include <syslog.h>
#endif
#include "GMS_malloc.h"
#include "GMS_indices.h"
#include "GMS_common.h"
#include "GMS_error_macros.h"

gms::math::
GrassScattererAVX::GrassScattererAVX() {

     m_gsc.nplants    = -1;
     m_gsc.nsteps     = -1;
     m_gsc.ordinal    = -1;
     m_gsc.param_npts = -1;
     m_gsc.tot_area   = 0.0f;
     m_gsc.lat        = 367.0f;
     m_gsc.lon        = 367.0f;
     m_gsc.elev       = -1.0f;
     m_gsc.Tah        = 0.0f;
     m_gsc.Tav        = 0.0f;
     m_gsc.epsilon    = {0.0f,0.0f};
     m_gsc.A          = NULL;
     m_gsc.moistness  = NULL;
     //
     m_gsc.xparam     = NULL;
     m_gsc.yparam     = NULL;
     m_gsc.zparam     = NULL;
     m_gsh.Polv[96]   = {};
     m_gsh.Polh[96]   = {};
     m_gsh.xang       = NULL;
     m_gsh.sin_xang   = NULL;
     m_gsh.cos_xang   = NULL;
     m_gsh.yang       = NULL;
     m_gsh.sin_yang   = NULL;
     m_gsh.cos_yang   = NULL;
}

gms::math::GrassScattererAVX::
GrassScattererAVX(const int32_t nsteps,
                  const int32_t ordinal,
		  const int32_t param_npts,
		  const float lat,
		  const float lon,
		  const float elev,
		  const std::complex<float> cepsilon) {

  using namespace gms::common;
  //generate random number of grass cylinders
  // per unit area.
  // Min value is 50 max value is 500
  const uint32_t lo   = 50U;
  const uint32_t hi   = 500U;
  int32_t result      = 0;
  uint32_t rand       = 0U;
  result = _rdrand32_step(&rand);
  if(!result) {
     	ABORT_ON_ERROR("GrassScattererAVX::GrassScattererAVX -- !!! _rdrand32_step failure !!! ", result)
  }
  if(rand < lo) rand = 50U;
  if(rand > hi) rand = 500U;
  m_gsc.nplants = static_cast<int32_t>(rand);
  m_gsc.nsteps = nsteps;
  m_gsc.ordinal = ordinal;
  m_gsc.param_npts = param_npts;
  m_gsc.tot_area  = 0.0f;
  m_gsc.lat       = lat;
  m_gsc.lon       = lon;
  m_gsc.elev      = elev;
  m_gsc.Tah       = 0.0f;
  m_gsc.Tav       = 0.0f;
  m_gsc.epsilon   = cepsilon;
  m_gsc.A         = gms_efmalloca(static_cast<size_t>(m_gsc.nplants),64);
  m_gsc.moistness = gms_eimalloca4(static_cast<size_t>(m_gsc.nplants),64);
  m_gsc.xparam    = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.param_npts),64);
  m_gsc.yparam    = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.param_npts),64);
  m_gsc.zparam    = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.param_npts),64);
  m_gsh.Polv[96]  = {0.0f,0.0f};
  m_gsh.Polh[96]  = {0.0f,0.0f};
  m_gsh.xang      = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.nsteps*m_gsc.nplants),64);
  m_gsh.sin_xang  = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.nsteps*m_gsc.nplants),64);
  m_gsh.cos_xang  = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.nsteps*m_gsc.nplants),64);
  m_gsh.yang      = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.nsteps*m_gsc.nplants),64);
  m_gsh.sin_yang  = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.nsteps*m_gsc.nplants),64);
  m_gsh.cos_yang  = gms_avxvec8_emalloca(static_cast<size_t>(m_gsc.nsteps*m_gsc.nplants),64);
}

gms::math::GrassScattererAVX
::~GrassScattererAVX() {

  _mm_free(m_gsc.A);
  m_gsc.A = NULL;
  _mm_free(m_gsc.moistness);
  m_gsc.moistness = NULL;
  _mm_free(m_gsc.xparam);
  m_gsc.xparam = NULL;
  _mm_free(m_gsc.yparam);
  m_gsc.yparam = NULL;
  _mm_free(m_gsc.zparam);
  m_gsc.zparam = NULL;
  _mm_free(m_gsh.xang);
  m_gsh.xang = NULL;
  _mm_free(m_gsh.sin_xang);
  m_gsh.sin_xang = NULL;
  _mm_free(m_gsh.cos_xang);
  m_gsh.cos_xang = NULL;
  _mm_free(m_gsh.yang);
  m_gsh.yang = NULL;
  _mm_free(m_gsh.sin_yang);
  m_gsh.sin_yang = NULL;
  _mm_free(m_gsh.cos_yang);
  m_gsh.cos_yang = NULL;
}

void
gms::math::GrassScattererAVX::
SetGrassMoistnessMask() {

  std::random_device rd;
  // Memory first touch
  for(int32_t i = 0; i != m_gsc.nplants; ++i) {
      m_gsc.moistness[i] = 0U;
  }
  std::mt19937 rgen(rd());
  std::uniform_int_distribution<> distro(0,1);
  for(int32_t i = 0; i != m_gsc.nplants; ++i) {
      m_gsc.moistness[i] = distro(rgen);
  }
}

void
gms::math::GrassScattererAVX::
ComputeGrassParamEq_ymm8r4(const char * __restrict pmc_event1,
			   const char * __restrict pmc_event2,
			   const char * __restrict pmc_event3,
			   const char * __restrict pmc_event4) {
                           
     struct _T0_ {
       AVXVec8 vtheta0;
       AVXVec8 vtheta1;
       AVXVec8 vtheta2;
       AVXVec8 vtheta3;
     } __ATTR_ALIGN__(64) t0;

     struct _T1_ {
       AVXVec8 vthinc0;
       AVXVec8 vthinc1;
       AVXVec8 vthinc2;
       AVXVec8 vthinc3;
     } __ATTR_ALIGN__(64) t1;

     struct _T2_ {
       AVXVec8 vhinc0;
       AVXVec8 vhinc1;
       AVXVec8 vhinc2;
       AVXVec8 vhinc3;
     } __ATTR_ALIGN__(64) t2;

     struct _T3_ {
       AVXVec8 vhinit0;
       AVXVec8 vhinit1;
       AVXVec8 vhinit2;
       AVXVec8 vhinit3;
     } __ATTR_ALIGN__(64) t3;

     struct _T4_ {
       AVXVec8 vNPTS;
       //
       AVXVec8 tmp1;
       AVXVec8 tmp2;
     } __ATTR_ALIGN__(64) t4;

     struct _T5_ {
       AVXVec8 tvrad;
       AVXVec8 tvz;
     } __ATTR_ALIGN__(64) t5;
     AVXVec8 rScale = AVXVec8{10.0f}; // unit of mm
     AVXVec8 hScale = AVXVec8{100.0f}; // unit of cm
     const int64_t totpts = static_cast<int64_t>(m_gsc.nplants*m_gsc.param_npts);
     const static float n2PI = 6.283185307179586f;
     std::clock_t seedr,seedz;
       //Locals first-touch
     t0.vtheta0 = ZERO;
     t0.vtheta1 = ZERO;
     t0.vtheta2 = ZERO;
     t0.vtheta3 = ZERO;
     t1.vthinc0 = ZERO;
     t1.vthinc1 = ZERO;
     t1.vthinc2 = ZERO;
     t1.vthinc3 = ZERO;
     t2.vhinc0  = ZERO;
     t2.vhinc1  = ZERO;
     t2.vhinc2  = ZERO;
     t2.vhinc3  = ZERO;
     t3.vhinit0 = ZERO;
     t3.vhinit1 = ZERO;
     t3.vhinit2 = ZERO;
     t3.vhinit3 = ZERO;
     t4.vNPTS   = ZERO;
     
     t4.tmp1    = ZERO;
     t4.tmp2    = ZERO;
     t4.vNPTS   = AVXVec8{static_cast<float>(m_tsc.param_npts)};
     
     t4.tmp1    = PI/t4.vNPTS;
     t5.tvrad   = ZERO;
     t5.tvz     = ZERO;
     // Memory first touch
     gms::common::avxvec8_init_unroll8x(&m_gsc.xparam[0],
                                        totpts,
					ZERO);
     t1.vthinc0 = t4.tmp1;
     t1.vthinc0 += VINC0;
     t1.vthinc1 = t4.tmp1;
     t1.vthinc1 += VINC1;
     gms::common::avxvec8_init_unroll8x(&m_gsc.yparam[0],
                                        totpts,
					ZERO);
     t1.vthinc2 = t4.tmp1;
     t1.vthinc2 += VINC2;
     t1.vthinc3 = t4.tmp1;
     t1.vthinc3 += VINC3;
     gms::common::avxvec8_init_unroll8x(&m_gsc.zparam[0],
                                        totpts,
					ZERO);
     // Loop over grass cylinders
#if defined __ICC || defined __INTEL_COMPILER
    
     __assume_aligned(m_gsc.A,64);
#elif defined __GNUC__ && !defined __INTEL_COMPILER
    
     m_gsc.A = (float*)__builtin_assume_aligned(m_gsc.A,64);
#endif
#if (SAMPLE_HW_PMC) == 1
          
            
	      // For now -- only single batch of 4 events is supported
	      const PFC_CNT ZERO_CNT[7] = {0,0,0,0,0,0,0};
	      PFC_CNT CNT[7] = {0,0,0,0,0,0,0};
	      PFC_CFG CFG[7] = {2,2,2,0,0,0,0};
	      CFG[3] = pfcParseCfg(pmc_event1);
	      CFG[4] = pfcParseCfg(pmc_event2);
	      CFG[5] = pfcParseCfg(pmc_event3);
	      CFG[6] = pfcParseCfg(pmc_event4);
	      // Reconfigure PMC and clear their count
	      pfcWrCfgs(0,7,CFG);
	      pfcWrCnts(0,7,ZERO_CNT);
	      memset(CNT,0,sizeof(CNT));
	      // Hot section
	      PFCSTART(CNT);
        
#endif
     for(int32_t i = 0; i != m_gsc.nplants; ++i) {
         seedr = std::clock();
	 auto rand_r = std::bind(std::uniform_real_distribution<float>(0.5f,0.9f),
	                         std::mt19937(seedr));
	 const float rtemp = rand_r();
	 //
	 t5.tvrad = rScale * rtemp;
	 seedz = std::clock();
	 auto rand_z = std::bind(std::uniform_real_distribution<float>(0.3f,0.5f),
	                         std::mt19937(seedz));
	 const float ztemp = rand_z();
	 const float c0 = n2PI*rtemp;
	 //
	 // Compute cylinders surface.
	 m_gsc.A[i] = c0*ztemp+c0*rtemp;
	 t5.tvz   = hScale * ztemp;
	 t4.tmp2  = t4.vNPTS;
	 t2.vhinc0 = t4.tmp2;
	 t2.vhinc0 += VINC0;
	 t2.vhinc1 = t4.tmp2;
	 t2.vhinc1 += VINC1;
	 t2.vhinc2 = t4.tmp2;
	 t2.vhinc2 += VINC2;
	 t2.vhinc3 = t4.tmp2;
	 t2.vhinc3 += VINC3;
	 t0.vtheta0 = ZERO;
	 t3.vhinit0 = ZERO;
	 t0.vtheta1 = ZERO;
	 t3.vhinit1 = ZERO;
	 t0.vtheta2 = ZERO;
	 t3.vhinit2 = ZERO;
	 t0.vtheta3 = ZERO;
	 t3.vhinit3 = ZERO;
#if defined __ICC || defined __INTEL_COMPILER
	 __assume_aligned(m_gsc.xparam,64);
	 __assume_aligned(m_gsc.yparam,64);
	 __assume_aligned(m_gsc.zparam,64);
#pragma vector always
#pragma vector vectorlength(8)
#pragma code_align(32)
#elif defined __GNUC__ && !defined __INTEL_COMPILER
        m_gsc.xparam = (AVXVec8*)__builtin_assume_aligned(m_gsc.xparam,64);
	m_gsc.yparam = (AVXVec8*)__builtin_assume_aligned(m_gsc.yparam,64);
	m_gsc.zparam = (AVXVec8*)__builtin_assume_aligned(m_gsc.zparam,64);
#pragma omp simd aligned(m_gsc.xparam,m_gsc.yparam,m_gsc.zparam:64)
#endif
        for(int32_t j = 0; j != m_gsc.param_npts-3; j += 4) {
             // for every point -- do...
	    t0.vtheta0 = t0.vtheta0+t1.vthinc0;
	    m_gsc.xparam[Ix2D(i,m_gsc.param_npts,j+0)] = t5.tvrad*cos(t0.vtheta0);
	    m_gsc.yparam[Ix2D(i,m_gsc.param_npts,j+0)] = t5.tvrad*sin(t0.vtheta0);
	    t3.vhinit0 = t3.vhinit0+t2.vhinc0;
	    m_gsc.zparam[Ix2D(i,m_gsc.param_npts,j+0)] = t3.vhinit0;
	    t0.vtheta1 = t0.vtheta1+t1.vthinc1;
	    m_gsc.xparam[Ix2D(i,m_gsc.param_npts,j+1)] = t5.tvrad*cos(t0.vtheta1);
	    m_gsc.yparam[Ix2D(i,m_gsc.param_npts,j+1)] = t5.tvrad*sin(t0.vtheta1);
	    t3.vhinit1 = t3.vhinit1+t2.vhinc1;
	    m_gsc.zparam[Ix2D(i,m_gsc.param_npts,j+1)] = t3.vhinit1;
	    t0.vtheta2 = t0.vtheta2+t1.vthinc2;
	    m_gsc.xparam[Ix2D(i,m_gsc.param_npts,j+2)] = t5.tvrad*cos(t0.vtheta2);
	    m_gsc.yparam[Ix2D(i,m_gsc.param_npts,j+2)] = t5.tvrad*sin(t0.vtheta2);
	    t3.vhinit2 = t3.vhinit2+t2.vhinc2;
	    m_gsc.zparam[Ix2D(i,m_gsc.param_npts,j+2)] = t3.vhinit2;
	    t0.vtheta3 = t0.vtheta3+t1.vthinc3;
	    m_gsc.xparam[Ix2D(i,m_gsc.param_npts,j+3)] = t5.tvrad*cos(t0.vtheta3);
	    m_gsc.yparam[Ix2D(i,m_gsc.param_npts,j+3)] = t5.tvrad*sin(t0.vtheta3);
	    t3.vhinit3 = t3.vhinit3+t2.vhinc3;
	    m_gsc.zparam[Ix2D(i,m_gsc.param_npts,j+3)] = t3.vhinit3;
	}
     }
#if (SAMPLE_HW_PMC) == 1
            PFCEND(CNT);
	    pfcRemoveBias(CNT,1);
	  
	    syslog(LOG_INFO,"%-10s:\n", __PRETTY_FUNCTION__);
	    syslog(LOG_INFO, "*************** Hardware Counters -- Dump Begin **************");
	    syslog(LOG_INFO,"Instructions Issued                  : %20lld\n", (signed long long)CNT[0]);
	    syslog(LOG_INFO,"Unhalted core cycles                 : %20lld\n", (signed long long)CNT[1]);
	    syslog(LOG_INFO,"Unhalted reference cycles            : %20lld\n", (signed long long)CNT[2]);
	    syslog(LOG_INFO,"%-37s: %20lld\n", pmc_event1                    , (signed long long)CNT[3]);
	    syslog(LOG_INFO,"%-37s: %20lld\n", pmc_event2                    , (signed long long)CNT[4]);
	    syslog(LOG_INFO,"%-37s: %20lld\n", pmc_event3                    , (signed long long)CNT[5]);
	    syslog(LOG_INFO,"%-37s: %20lld\n", pmc_event4                    , (signed long long)CNT[6]);
	    syslog(LOG_INFO, "*************** Hardware Counters -- Dump End   **************");
#endif	
}

#include <math.h>
void
gms::math::GrassScattererAVX::
ComputeGrassHVPolarization(
			   const float gamma,
			   const float ah,   
			   const float av) {
     const static float n2PI = 6.283185307179586f;
     const static float n1rad = 0.0174444444444444444f;
     std::complex<float> dielectric = {0.0f,0.0f};
     std::complex<float> ch0        = {0.0f,0.0f};
     //
     std::complex<float> cv0        = {0.0f,0.0f};
     //
     const float n28PI = 14.0f*TWOPI;
     const float K = n2PI/gamma;
     const float sarea = 0.0f;
     float term1 = 0.0f;
     float term2 = 0.0f;
     float term3 = 0.0;
     float term4 = 0.0f;
     float term5 = 0.0f;
     float t     = 0.0f;
     float re    = 0.0f;
     float im    = 0.0f;
     float c0    = 0.0f;
     float c1    = 0.0f;
#if defined __ICC || defined __INTEL_COMPILER
     __assume_aligned(m_gsc.A,64);
#elif defined __GNUC__ && !defined __INTEL_COMPILER
     m_gsc.A = (float*)__builtin_assume_aligned(m_gsc.A,64);
#endif
#if defined __ICC || defined __INTEL_COMPILER
#pragma code_align(32)

#endif
     for(int32_t i = 0; i != m_gsc.nplants-3; i += 4) {
     
         float t0 += m_gsc.A[i+0];
	 float t1 += m_gsc.A[i+1];
	 float t2 += m_gsc.A[i+2];
	 float t3 += m_gsc.A[i+3]
	// const float tot = t0+t1;
	// m_gsc.tot_area += tot;
     }
     m_gsc.tot_area = t0+t1+t2+t3;
     sarea = m_gsc.tot_area*m_gsc.tot_area;
     term1 = static_cast<float>(m_gsc.nplants)*asqrt*K*K;
     c0 = m_gsc.epsilon.real();
     t = 1.0f/(1.0f+c0);
     c1 = m_gsc.epsilon.imag();
     re = c0-1.0f*c0-1.0f;
     im = c1*c1;
     dielectric = {c0,c1};
     term2 = 3.0f+16.0f*t+96.0f*t*t;
     term3 = 3.0f*((ah/K)*(ah/K));
     term4 = 12+8.0f*t-64.0f*t*t;
     term5 = 3.0f*((av/K)*(av/K));
     for(int32_t i = 0; i != 90; ++i) {
         // Full grazing angle sweep.
	 const float theta = n1rad*static_cast<float>(i);
	 const float t0    = sin(theta);
	 const float a0    = 4.0f*(1.0f+2.0f*t0*t0);
	 ch0  = dielectric*term2;
	 const float t1    = term1/(n28PI*t0);
	 const float t2    = term3+a0;
	 m_gsh.Polh[i]     = t1*ch0/t2;
	 const float t3    = cos(theta);
	 const float t4    = term2+t3*t3*term4;
	 cv1 = dielectric*t4;
	 const float t5    = term5+a0;
	 m_gsh.Polv[i]     = t1*cv1/t5;
     }
     
}



