#include "SamplingStrategy.h"
#include "Ray.h"
#include "Intersectable.h"
#include "Preprocessor.h"
#include "ViewCellsManager.h"
#include "AxisAlignedBox3.h"
#include "RssTree.h"
#include "Vector2.h"
#include "RndGauss.h"
#include "Mutation.h"

#ifdef GTP_INTERNAL
#include "ArchModeler2MLRT.hxx"
#endif

namespace GtpVisibilityPreprocessor {

#define MUTATION_USE_CDF 0
#define USE_SILHOUETTE_MUTATIONS 0
#define EVALUATE_MUTATION_STATS 1

void
MutationBasedDistribution::Update(VssRayContainer &vssRays)
{
  //  for (int i=0; i < mRays.size(); i++)
  //	cout<<mRays[i].mMutations<<" ";
  //  cout<<endl;
  cerr<<"Muattion update..."<<endl;
  cerr<<"rays = "<<mRays.size()<<endl;
  if (mRays.size()) {
	cerr<<"Oversampling factors = "<<
	  GetEntry(0).mMutations<<" "<<
	  GetEntry(1).mMutations<<" "<<
	  GetEntry(2).mMutations<<" "<<
	  GetEntry(3).mMutations<<" "<<
	  GetEntry(4).mMutations<<" "<<
	  GetEntry(5).mMutations<<" ... "<<
	  GetEntry(mRays.size()-6).mMutations<<" "<<
	  GetEntry(mRays.size()-5).mMutations<<" "<<
	  GetEntry(mRays.size()-4).mMutations<<" "<<
	  GetEntry(mRays.size()-3).mMutations<<" "<<
	  GetEntry(mRays.size()-2).mMutations<<" "<<
	  GetEntry(mRays.size()-1).mMutations<<endl;
  }
  int contributingRays = 0;

  int mutationRays = 0;
  int dummyNcMutations = 0;
  int dummyCMutations = 0;
  
  for (int i=0; i < vssRays.size(); i++) {
	if (vssRays[i]->mPvsContribution) {
	  // reset the counter of unsuccsseful mutation for a generating ray (if it exists)
	  if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
		  vssRays[i]->mGeneratorId != -1
		  ) {
		mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
#if EVALUATE_MUTATION_STATS 
		mutationRays++;
		
		Intersectable *newObject =
		  mPreprocessor.mViewCellsManager->GetIntersectable(
															*vssRays[i],
															true);

		Intersectable *oldObject =
		  mPreprocessor.mViewCellsManager->GetIntersectable(
															*mRays[vssRays[i]->mGeneratorId].mRay,
															true);
		
		if (oldObject == newObject)
		  dummyCMutations++;
#endif
	  }
	  contributingRays++;
	  if (mRays.size() < mMaxRays) {
		VssRay *newRay = new VssRay(*vssRays[i]);
		// add this ray
		newRay->Ref();
		mRays.push_back(RayEntry(newRay));
	  } else {
		// unref the old ray
		*mRays[mBufferStart].mRay = *vssRays[i];
		mRays[mBufferStart].mMutations = 0;
		//		mRays[mBufferStart] = RayEntry(newRay);
		mBufferStart++;
		if (mBufferStart >= mMaxRays)
		  mBufferStart = 0;
	  }
	} else {
#if EVALUATE_MUTATION_STATS 
	  if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
		  vssRays[i]->mGeneratorId != -1
		  ) {
		mutationRays++;

		Intersectable *newObject =
		  mPreprocessor.mViewCellsManager->GetIntersectable(
															*vssRays[i],
															true);

		Intersectable *oldObject =
		  mPreprocessor.mViewCellsManager->GetIntersectable(
															*mRays[vssRays[i]->mGeneratorId].mRay,
															true);

		if (oldObject == newObject)
		  dummyNcMutations++;
	  }
#endif
	}
  }

  if (mutationRays) {
	cout<<"Mutated rays:"<<mutationRays<<endl;
	cout<<"Dummy mutations ratio:"<<100.0f*(dummyCMutations + dummyNcMutations)/
	  (float)mutationRays<<"%"<<endl;
	cout<<"Dummy NC mutations ratio:"<<100.0f*dummyNcMutations/(float)mutationRays<<"%"<<endl;
	cout<<"Dummy C mutations ratio:"<<100.0f*dummyCMutations/(float)mutationRays<<"%"<<endl;
  }

  float pContributingRays = contributingRays/(float)vssRays.size();
  float importance = 1.0f/(pContributingRays + 1e-5);
  // set this values for last contributingRays
  int index = mBufferStart - 1;
  
  for (int i=0; i < contributingRays; i++, index--) {
	if (index < 0)
	  index = mRays.size()-1;
	mRays[index].mImportance = importance;
  }

#if MUTATION_USE_CDF
  // compute cdf
  mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mMutations+1);
  for (int i=1; i < mRays.size(); i++) 
	mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mMutations+1);
  
  float scale = 1.0f/mRays[i-1].mCdf;
  for (i=0; i < mRays.size(); i++) {
	mRays[i].mCdf *= scale;
  }
#endif
  
  cout<<"Importance = "<<
	GetEntry(0).mImportance<<" "<<
	GetEntry(mRays.size()-1).mImportance<<endl;
  
  cerr<<"Mutation update done."<<endl;
}


Vector3
MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray,
												 const Vector3 &U,
												 const Vector3 &V,
												 const Vector2 vr2,
												 const float radius
												 )
{
#if 0
  Vector3 v;
  if (d.DrivingAxis() == 0)
	v = Vector3(0, r[0]-0.5f, r[1]-0.5f);
  else
	if (d.DrivingAxis() == 1)
	  v = Vector3(r[0]-0.5f, 0, r[1]-0.5f);
	else
	  v = Vector3(r[0]-0.5f, r[1]-0.5f, 0);
  return v*(2*radius);
#endif
#if 0
  return (U*(r[0] - 0.5f) + V*(r[1] - 0.5f))*(2*radius);
#endif


  // Output random variable
  Vector2 gaussvec2; 
  
  // Here we apply transform to gaussian, so 2D bivariate
  // normal distribution
  //  float sigma = ComputeSigmaFromRadius(radius);
  float sigma = radius;
  GaussianOn2D(vr2,
			   sigma, // input
			   gaussvec2);  // output

	
    // Here we tranform the point correctly to 3D space using base
    // vectors of the 3D space defined by the direction
    Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
	
	//	cout<<shift<<endl;
	return shift;
}

Vector3
MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray,
													  const Vector3 &U,
													  const Vector3 &V,
													  const Vector2 vr2,
													  const float radius
													  )
{
#if 0
  Vector3 v;
  // mutate the termination
  if (d.DrivingAxis() == 0)
	v = Vector3(0, r[2]-0.5f, r[3]-0.5f);
  else
	if (d.DrivingAxis() == 1)
	  v = Vector3(r[2]-0.5f, 0, r[3]-0.5f);
	else
	  v = Vector3(r[2]-0.5f, r[3]-0.5f, 0);
  
  //   Vector3 nv;
  
  //   if (Magnitude(v) > Limits::Small)
  // 	nv = Normalize(v);
  //   else
  // 	nv = v;
  
  //  v = nv*size + v*size;

  return v*(4.0f*radius);
#endif
#if 0
  return (U*(vr2.xx - 0.5f) + V*(vr2.yy - 0.5f))*(4.0f*radius);
#endif
  Vector2 gaussvec2; 
#if 1
  float sigma = radius;
  GaussianOn2D(vr2,
			   sigma, // input
			   gaussvec2);  // output
  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
  //	cout<<shift<<endl;
  return shift;
#endif
#if 0
  // Here we estimate standard deviation (sigma) from radius
  float sigma = 1.1f*ComputeSigmaFromRadius(radius);
  Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1));
  PolarGaussianOnDisk(vr3,
					  sigma,
					  radius, // input
					  gaussvec2); // output
  
  // Here we tranform the point correctly to 3D space using base
  // vectors of the 3D space defined by the direction
  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
  
  //	cout<<shift<<endl;
  return shift;
#endif
}




bool
MutationBasedDistribution::GenerateSample(SimpleRay &sray)
{

  if (mRays.size() == 0) {
	float rr[5];
	// use direction based distribution
	Vector3 origin, direction;
	static HaltonSequence halton;
	
	halton.GetNext(5, rr);
	mPreprocessor.mViewCellsManager->GetViewPoint(origin,
												  Vector3(rr[0], rr[1], rr[2]));
	

	direction = UniformRandomVector(rr[3], rr[4]);
	
	const float pdf = 1.0f;
	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
	sray.mGeneratorId = -1;

	return true;
  } 

  int index;

#if !MUTATION_USE_CDF
  // get tail of the buffer
  index = (mLastIndex+1)%mRays.size();
  if (mRays[index].GetSamplingFactor() >
	  mRays[mLastIndex].GetSamplingFactor()) {
	// search back for index where this is valid
	index = (mLastIndex - 1 + mRays.size())%mRays.size();
	for (int i=0; i < mRays.size(); i++) {
	  
	  //	  if (mRays[index].mMutations > mRays[mLastIndex].mMutations)
	  //		break;
	  if (mRays[index].GetSamplingFactor() >
		  mRays[mLastIndex].GetSamplingFactor() )
		break;
	  index = (index - 1 + mRays.size())%mRays.size();
	}
	// go one step back
	index = (index+1)%mRays.size();
  }
#else
  static HaltonSequence iHalton;
  iHalton.GetNext(1, rr);
  //rr[0] = RandomValue(0,1);
  // use binary search to find index with this cdf
  int l=0, r=mRays.size()-1;
  while(l<r) {
	int i = (l+r)/2;
	if (rr[0] < mRays[i].mCdf )
	  r = i;
	else
	  l = i+1;
  }
  index = l;
  //  if (rr[0] >= mRays[r].mCdf)
  //	index = r;
  //  else
  //	index = l;
	
	
#endif
  //  cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl;

  mLastIndex = index;

#if USE_SILHOUETTE_MUTATIONS
  return GenerateSilhouetteMutation(index, sray);
#else  
  return GenerateMutation(index, sray);
#endif
  
}





bool
MutationBasedDistribution::GenerateMutationCandidate(const int index,
													 SimpleRay &sray,
													 Intersectable *object,
													 const AxisAlignedBox3 &box
													 )
{
  float rr[4];

  VssRay *ray = mRays[index].mRay;

  mRays[index].mHalton.GetNext(4, rr);
  
  // mutate the origin
  Vector3 d = ray->GetDir();
  
  float objectRadius = 0.5f*Magnitude(box.Diagonal());
  //  cout<<objectRadius<<endl;
  if (objectRadius < Limits::Small)
	return false;
  
  // Compute right handed coordinate system from direction
  Vector3 U, V;
  Vector3 nd = Normalize(d);
  nd.RightHandedBase(U, V);
  
  Vector3 origin = ray->mOrigin;
  Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center();
  
  float radiusExtension = 0.05f;
  //  + mRays[index].mMutations/50.0f;
  
  origin += ComputeOriginMutation(*ray, U, V,
								  Vector2(rr[0], rr[1]),
								  objectRadius*radiusExtension);
  
  termination += ComputeTerminationMutation(*ray, U, V,
											Vector2(rr[2], rr[3]),
											objectRadius*radiusExtension);
  
  Vector3 direction = termination - origin;
  
  if (Magnitude(direction) < Limits::Small)
	return false;
  
  // shift the origin a little bit
  origin += direction*0.5f;
  
  direction.Normalize();
  
  // $$ jb the pdf is yet not correct for all sampling methods!
  const float pdf = 1.0f;
  
  sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
  sray.mGeneratorId = index;
}

bool
MutationBasedDistribution::GenerateMutation(const int index, SimpleRay &sray)
{
  VssRay *ray = mRays[index].mRay;

  Intersectable *object = mPreprocessor.mViewCellsManager->GetIntersectable(
																			*ray,
																			true);
  
  AxisAlignedBox3 box = object->GetBox();
  
  if (GenerateMutationCandidate(index, sray, object, box)) {
    mRays[index].mMutations++;
	return true;
  }
  return false;
}

bool
MutationBasedDistribution::GenerateSilhouetteMutation(const int index, SimpleRay &sray)
{
#ifndef GTP_INTERNAL
  return GenerateMutation(index, sray);
#else
  const int packetSize = 4;
  const int maxTries = 8;
  
  static int hit_triangles[16];
  static float dist[16];
  
  SimpleRay mutationCandidates[packetSize];
  int candidates = 0;

  VssRay *ray = mRays[index].mRay;
  
  Intersectable *object = mPreprocessor.mViewCellsManager->GetIntersectable(
																			*ray,
																			true);
  
  AxisAlignedBox3 box = object->GetBox();

  int id = 0;
  int silhouetteRays = 0;
  int tries = 0;
  while (silhouetteRays == 0 && tries < maxTries) {
	for (candidates = 0; candidates < packetSize && tries < maxTries; tries++)
	  if (GenerateMutationCandidate(index, mutationCandidates[candidates], object, box))
		candidates++;

	if (candidates < packetSize)
	  break;
	
	//	cout<<candidates<<endl;
	// cast rays to find silhouette edge
	for (int i=0; i < packetSize; i++)
	  mlrtaStoreRayAS4(&mutationCandidates[i].mOrigin.x,
					   &mutationCandidates[i].mDirection.x,
					   i);

	mlrtaTraverseGroupAS4(&box.Min().x,
						  &box.Max().x,
						  hit_triangles,
						  dist);
	
	for (int i=0; i < packetSize; i++)
	  if (hit_triangles[i] == -1) {
		silhouetteRays++;
		id = i;
		break;
	  }
  }

  if (candidates == 0)
	return false;
  
  //  cout<<id<<endl;
  // cout<<tries<<endl;
  sray = mutationCandidates[id];
  mRays[index].mMutations++;

  return true;
#endif
}

  


MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor
													 ) :
  SamplingStrategy(preprocessor)
{
  mType = MUTATION_BASED_DISTRIBUTION;
  mBufferStart = 0;
  mMaxRays = 500000;
  mRays.reserve(mMaxRays);
  mOriginMutationSize = 10.0f;
  mLastIndex = 0;
  //  mOriginMutationSize = Magnitude(preprocessor.mViewCellsManager->
  //								  GetViewSpaceBox().Diagonal())*1e-3;
  
}


}