// Copyright 2009-2021 Intel Corporation
// SPDX-License-Identifier: Apache-2.0

#pragma once

#include "../builders/primrefgen.h"
#include "../builders/heuristic_spatial.h"
#include "../builders/splitter.h"

#include "../../common/algorithms/parallel_for_for.h"
#include "../../common/algorithms/parallel_for_for_prefix_sum.h"

#define DBG_PRESPLIT(x)   
#define CHECK_PRESPLIT(x) 

#define GRID_SIZE 1024
#define MAX_PRESPLITS_PER_PRIMITIVE_LOG 5
#define MAX_PRESPLITS_PER_PRIMITIVE (1<<MAX_PRESPLITS_PER_PRIMITIVE_LOG)
#define PRIORITY_CUTOFF_THRESHOLD 1.0f
#define PRIORITY_SPLIT_POS_WEIGHT 1.5f

namespace embree
{  
  namespace isa
  {

    struct PresplitItem
    {
      union {
        float priority;    
        unsigned int data;
      };
      unsigned int index;
      
      __forceinline operator unsigned() const
      {
	return reinterpret_cast<const unsigned&>(priority);
      }
      __forceinline bool operator < (const PresplitItem& item) const
      {
	return (priority < item.priority);
      }

      template<typename Mesh>
      __forceinline static float compute_priority(const PrimRef &ref, Scene *scene, const Vec2i &mc)
      {
	const unsigned int geomID = ref.geomID();
	const unsigned int primID = ref.primID();
	const float area_aabb  = area(ref.bounds());
	const float area_prim  = ((Mesh*)scene->get(geomID))->projectedPrimitiveArea(primID);
        const unsigned int diff = 31 - lzcnt(mc.x^mc.y);
        assert(area_prim <= area_aabb);
        //const float priority = powf((area_aabb - area_prim) * powf(PRIORITY_SPLIT_POS_WEIGHT,(float)diff),1.0f/4.0f);   
        const float priority = sqrtf(sqrtf( (area_aabb - area_prim) * powf(PRIORITY_SPLIT_POS_WEIGHT,(float)diff) ));
        assert(priority >= 0.0f && priority < FLT_LARGE);
	return priority;      
      }

    
    };

    inline std::ostream &operator<<(std::ostream &cout, const PresplitItem& item) {
      return cout << "index " << item.index << " priority " << item.priority;    
    };

    template<typename SplitterFactory>    
      void splitPrimitive(SplitterFactory &Splitter,
                          const PrimRef &prim,
                          const unsigned int geomID,
                          const unsigned int primID,
                          const unsigned int split_level,
                          const Vec3fa &grid_base, 
                          const float grid_scale,
                          const float grid_extend,
                          PrimRef subPrims[MAX_PRESPLITS_PER_PRIMITIVE],
                          unsigned int& numSubPrims)
    {
      assert(split_level <= MAX_PRESPLITS_PER_PRIMITIVE_LOG);
      if (split_level == 0)
      {
        assert(numSubPrims < MAX_PRESPLITS_PER_PRIMITIVE);
        subPrims[numSubPrims++] = prim;
      }
      else
      {
        const Vec3fa lower = prim.lower;
        const Vec3fa upper = prim.upper;
        const Vec3fa glower = (lower-grid_base)*Vec3fa(grid_scale)+Vec3fa(0.2f);
        const Vec3fa gupper = (upper-grid_base)*Vec3fa(grid_scale)-Vec3fa(0.2f);
        Vec3ia ilower(floor(glower));
        Vec3ia iupper(floor(gupper));

        /* this ignores dimensions that are empty */
        iupper = (Vec3ia)(select(vint4(glower) >= vint4(gupper),vint4(ilower),vint4(iupper)));

        /* compute a morton code for the lower and upper grid coordinates. */
        const unsigned int lower_code = bitInterleave(ilower.x,ilower.y,ilower.z);
        const unsigned int upper_code = bitInterleave(iupper.x,iupper.y,iupper.z);
			
        /* if all bits are equal then we cannot split */
        if(unlikely(lower_code == upper_code))
        {
          assert(numSubPrims < MAX_PRESPLITS_PER_PRIMITIVE);
          subPrims[numSubPrims++] = prim;
          return;
        }
		    
        /* compute octree level and dimension to perform the split in */
        const unsigned int diff = 31 - lzcnt(lower_code^upper_code);
        const unsigned int level = diff / 3;
        const unsigned int dim   = diff % 3;
      
        /* now we compute the grid position of the split */
        const unsigned int isplit = iupper[dim] & ~((1<<level)-1);
			    
        /* compute world space position of split */
        const float inv_grid_size = 1.0f / GRID_SIZE;
        const float fsplit = grid_base[dim] + isplit * inv_grid_size * grid_extend;

        assert(prim.lower[dim] <= fsplit &&
               prim.upper[dim] >= fsplit);
		
        /* split primitive */
        const auto splitter = Splitter(prim);
        BBox3fa left,right;
        splitter(prim.bounds(),dim,fsplit,left,right);
        assert(!left.empty());
        assert(!right.empty());

			    
        splitPrimitive(Splitter,PrimRef(left ,geomID,primID),geomID,primID,split_level-1,grid_base,grid_scale,grid_extend,subPrims,numSubPrims);
        splitPrimitive(Splitter,PrimRef(right,geomID,primID),geomID,primID,split_level-1,grid_base,grid_scale,grid_extend,subPrims,numSubPrims);
      }
    }
    
    
    template<typename Mesh, typename SplitterFactory>    
      PrimInfo createPrimRefArray_presplit(Geometry* geometry, unsigned int geomID, size_t numPrimRefs, mvector<PrimRef>& prims, BuildProgressMonitor& progressMonitor)
    {
      ParallelPrefixSumState<PrimInfo> pstate;
      
      /* first try */
      progressMonitor(0);
      PrimInfo pinfo = parallel_prefix_sum( pstate, size_t(0), geometry->size(), size_t(1024), PrimInfo(empty), [&](const range<size_t>& r, const PrimInfo& base) -> PrimInfo {
	  return geometry->createPrimRefArray(prims,r,r.begin(),geomID);
	}, [](const PrimInfo& a, const PrimInfo& b) -> PrimInfo { return PrimInfo::merge(a,b); });

      /* if we need to filter out geometry, run again */
      if (pinfo.size() != numPrimRefs)
	{
	  progressMonitor(0);
	  pinfo = parallel_prefix_sum( pstate, size_t(0), geometry->size(), size_t(1024), PrimInfo(empty), [&](const range<size_t>& r, const PrimInfo& base) -> PrimInfo {
	      return geometry->createPrimRefArray(prims,r,base.size(),geomID);
	    }, [](const PrimInfo& a, const PrimInfo& b) -> PrimInfo { return PrimInfo::merge(a,b); });
	}
      return pinfo;	
    }
    
    __forceinline Vec2i computeMC(const Vec3fa &grid_base, const float grid_scale, const PrimRef &ref)
    {
      const Vec3fa lower = ref.lower;
      const Vec3fa upper = ref.upper;
      const Vec3fa glower = (lower-grid_base)*Vec3fa(grid_scale)+Vec3fa(0.2f);
      const Vec3fa gupper = (upper-grid_base)*Vec3fa(grid_scale)-Vec3fa(0.2f);
      Vec3ia ilower(floor(glower));
      Vec3ia iupper(floor(gupper));
      
      /* this ignores dimensions that are empty */
      iupper = (Vec3ia)select(vint4(glower) >= vint4(gupper),vint4(ilower),vint4(iupper));

      /* compute a morton code for the lower and upper grid coordinates. */
      const unsigned int lower_code = bitInterleave(ilower.x,ilower.y,ilower.z);
      const unsigned int upper_code = bitInterleave(iupper.x,iupper.y,iupper.z);
      return Vec2i(lower_code,upper_code);
    }

    template<typename Mesh, typename SplitterFactory>    
      PrimInfo createPrimRefArray_presplit(Scene* scene, Geometry::GTypeMask types, bool mblur, size_t numPrimRefs, mvector<PrimRef>& prims, BuildProgressMonitor& progressMonitor)
    {	
      static const size_t MIN_STEP_SIZE = 128;

      ParallelForForPrefixSumState<PrimInfo> pstate;
      Scene::Iterator2 iter(scene,types,mblur);

      /* first try */
      progressMonitor(0);
      pstate.init(iter,size_t(1024));
      PrimInfo pinfo = parallel_for_for_prefix_sum0( pstate, iter, PrimInfo(empty), [&](Geometry* mesh, const range<size_t>& r, size_t k, size_t geomID) -> PrimInfo {
	  return mesh->createPrimRefArray(prims,r,k,(unsigned)geomID);
	}, [](const PrimInfo& a, const PrimInfo& b) -> PrimInfo { return PrimInfo::merge(a,b); });
      
      /* if we need to filter out geometry, run again */
      if (pinfo.size() != numPrimRefs)
	{
	  progressMonitor(0);
	  pinfo = parallel_for_for_prefix_sum1( pstate, iter, PrimInfo(empty), [&](Geometry* mesh, const range<size_t>& r, size_t k, size_t geomID, const PrimInfo& base) -> PrimInfo {
	      return mesh->createPrimRefArray(prims,r,base.size(),(unsigned)geomID);
	    }, [](const PrimInfo& a, const PrimInfo& b) -> PrimInfo { return PrimInfo::merge(a,b); });
	}

      /* use correct number of primitives */
      size_t numPrimitives = pinfo.size();
      const size_t alloc_numPrimitives = prims.size(); 
      const size_t numSplitPrimitivesBudget = alloc_numPrimitives - numPrimitives;

      /* set up primitive splitter */
      SplitterFactory Splitter(scene);


      DBG_PRESPLIT(
        const size_t org_numPrimitives = pinfo.size();
        PRINT(numPrimitives);		
        PRINT(alloc_numPrimitives);		
        PRINT(numSplitPrimitivesBudget);
        );

      /* allocate double buffer presplit items */
      const size_t presplit_allocation_size = sizeof(PresplitItem)*alloc_numPrimitives;
      PresplitItem *presplitItem     = (PresplitItem*)alignedMalloc(presplit_allocation_size,64);
      PresplitItem *tmp_presplitItem = (PresplitItem*)alignedMalloc(presplit_allocation_size,64);

      /* compute grid */
      const Vec3fa grid_base    = pinfo.geomBounds.lower;
      const Vec3fa grid_diag    = pinfo.geomBounds.size();
      const float grid_extend   = max(grid_diag.x,max(grid_diag.y,grid_diag.z));		
      const float grid_scale    = grid_extend == 0.0f ? 0.0f : GRID_SIZE / grid_extend;

      /* init presplit items and get total sum */
      const float psum = parallel_reduce( size_t(0), numPrimitives, size_t(MIN_STEP_SIZE), 0.0f, [&](const range<size_t>& r) -> float {
          float sum = 0.0f;
          for (size_t i=r.begin(); i<r.end(); i++)
          {		
            presplitItem[i].index = (unsigned int)i;
            const Vec2i mc = computeMC(grid_base,grid_scale,prims[i]);
            /* if all bits are equal then we cannot split */
            presplitItem[i].priority = (mc.x != mc.y) ? PresplitItem::compute_priority<Mesh>(prims[i],scene,mc) : 0.0f;    
            /* FIXME: sum undeterministic */
            sum += presplitItem[i].priority;
          }
          return sum;
        },[](const float& a, const float& b) -> float { return a+b; });

      /* compute number of splits per primitive */
      const float inv_psum = 1.0f / psum;
      parallel_for( size_t(0), numPrimitives, size_t(MIN_STEP_SIZE), [&](const range<size_t>& r) -> void {
          for (size_t i=r.begin(); i<r.end(); i++)
          {
            if (presplitItem[i].priority > 0.0f)
            {
              const float rel_p = (float)numSplitPrimitivesBudget * presplitItem[i].priority * inv_psum;
              if (rel_p >= PRIORITY_CUTOFF_THRESHOLD) // need at least a split budget that generates two sub-prims
              {
                presplitItem[i].priority = max(min(ceilf(logf(rel_p)/logf(2.0f)),(float)MAX_PRESPLITS_PER_PRIMITIVE_LOG),1.0f);
                //presplitItem[i].priority = min(floorf(logf(rel_p)/logf(2.0f)),(float)MAX_PRESPLITS_PER_PRIMITIVE_LOG);
                assert(presplitItem[i].priority >= 0.0f && presplitItem[i].priority <= (float)MAX_PRESPLITS_PER_PRIMITIVE_LOG);
              }
              else
                presplitItem[i].priority = 0.0f;
            }
          }
        });

      auto isLeft = [&] (const PresplitItem &ref) { return ref.priority < PRIORITY_CUTOFF_THRESHOLD; };        
      size_t center = parallel_partitioning(presplitItem,0,numPrimitives,isLeft,1024);

      /* anything to split ? */
      if (center < numPrimitives)
      {
        size_t numPrimitivesToSplit = numPrimitives - center;
        assert(presplitItem[center].priority >= 1.0f);

        /* sort presplit items in ascending order */
        radix_sort_u32(presplitItem + center,tmp_presplitItem + center,numPrimitivesToSplit,1024);

        CHECK_PRESPLIT(
          parallel_for( size_t(center+1), numPrimitives, size_t(MIN_STEP_SIZE), [&](const range<size_t>& r) -> void {
              for (size_t i=r.begin(); i<r.end(); i++)
                assert(presplitItem[i-1].priority <= presplitItem[i].priority);
            });
          );

        unsigned int* primOffset0 = (unsigned int*)tmp_presplitItem;
        unsigned int* primOffset1 = (unsigned int*)tmp_presplitItem + numPrimitivesToSplit;

        /* compute actual number of sub-primitives generated within the [center;numPrimitives-1] range */
        const size_t totalNumSubPrims = parallel_reduce( size_t(center), numPrimitives, size_t(MIN_STEP_SIZE), size_t(0), [&](const range<size_t>& t) -> size_t {
            size_t sum = 0;
            for (size_t i=t.begin(); i<t.end(); i++)
            {	
              PrimRef subPrims[MAX_PRESPLITS_PER_PRIMITIVE];	
              assert(presplitItem[i].priority >= 1.0f);
              const unsigned int  primrefID = presplitItem[i].index;	
              const float prio              = presplitItem[i].priority;
              const unsigned int   geomID   = prims[primrefID].geomID();
              const unsigned int   primID   = prims[primrefID].primID();
              const unsigned int split_levels = (unsigned int)prio;
              unsigned int numSubPrims = 0;
              splitPrimitive(Splitter,prims[primrefID],geomID,primID,split_levels,grid_base,grid_scale,grid_extend,subPrims,numSubPrims);
              assert(numSubPrims);
              numSubPrims--; // can reuse slot 
              sum+=numSubPrims;
              presplitItem[i].data = (numSubPrims << MAX_PRESPLITS_PER_PRIMITIVE_LOG) | split_levels;
              primOffset0[i-center] = numSubPrims;
            }
            return sum;
          },[](const size_t& a, const size_t& b) -> size_t { return a+b; });
        
        /* if we are over budget, need to shrink the range */
        if (totalNumSubPrims > numSplitPrimitivesBudget) 
        {
          size_t new_center = numPrimitives-1;
          size_t sum = 0;
          for (;new_center>=center;new_center--)
          {
            const unsigned int numSubPrims = presplitItem[new_center].data >> MAX_PRESPLITS_PER_PRIMITIVE_LOG;
            if (unlikely(sum + numSubPrims >= numSplitPrimitivesBudget)) break;
            sum += numSubPrims;
          }
          new_center++;

          primOffset0 += new_center - center;
          numPrimitivesToSplit -= new_center - center;
          center = new_center;
          assert(numPrimitivesToSplit == (numPrimitives - center));
        }

        /* parallel prefix sum to compute offsets for storing sub-primitives */
        const unsigned int offset = parallel_prefix_sum(primOffset0,primOffset1,numPrimitivesToSplit,(unsigned int)0,std::plus<unsigned int>());
        assert(numPrimitives+offset <= alloc_numPrimitives);

        /* iterate over range, and split primitives into sub primitives and append them to prims array */		    
        parallel_for( size_t(center), numPrimitives, size_t(MIN_STEP_SIZE), [&](const range<size_t>& rn) -> void {
            for (size_t j=rn.begin(); j<rn.end(); j++)		    
            {
              PrimRef subPrims[MAX_PRESPLITS_PER_PRIMITIVE];
              const unsigned int  primrefID = presplitItem[j].index;	
              const unsigned int   geomID   = prims[primrefID].geomID();
              const unsigned int   primID   = prims[primrefID].primID();
              const unsigned int split_levels = presplitItem[j].data & ((unsigned int)(1 << MAX_PRESPLITS_PER_PRIMITIVE_LOG)-1);

              assert(split_levels);
              assert(split_levels <= MAX_PRESPLITS_PER_PRIMITIVE_LOG);
              unsigned int numSubPrims = 0;
              splitPrimitive(Splitter,prims[primrefID],geomID,primID,split_levels,grid_base,grid_scale,grid_extend,subPrims,numSubPrims);
              const size_t newID = numPrimitives + primOffset1[j-center];              
              assert(newID+numSubPrims-1 <= alloc_numPrimitives);
              prims[primrefID] = subPrims[0];
              for (size_t i=1;i<numSubPrims;i++)
                prims[newID+i-1] = subPrims[i];
            }
          });

        numPrimitives += offset;
        DBG_PRESPLIT(
          PRINT(pinfo.size());
          PRINT(numPrimitives);
          PRINT((float)numPrimitives/org_numPrimitives));                
      }
                
      /* recompute centroid bounding boxes */
      pinfo = parallel_reduce(size_t(0),numPrimitives,size_t(MIN_STEP_SIZE),PrimInfo(empty),[&] (const range<size_t>& r) -> PrimInfo {
          PrimInfo p(empty);
          for (size_t j=r.begin(); j<r.end(); j++)
            p.add_center2(prims[j]);
          return p;
        }, [](const PrimInfo& a, const PrimInfo& b) -> PrimInfo { return PrimInfo::merge(a,b); });
  
      assert(pinfo.size() == numPrimitives);
      
      /* free double buffer presplit items */
      alignedFree(tmp_presplitItem);		
      alignedFree(presplitItem);
      return pinfo;	
    }
  }
}