package edu.mit.husain.DSMCVR3D;

import java.io.Serializable;
import static edu.mit.husain.DSMCVR3D.myUtils.*;
import static java.lang.Math.*;


public class KDTree implements Serializable{
	private static final long serialVersionUID = 1L;
	private static float BIG=Float.MAX_VALUE;
	private Box root;
	final double[] x,y,z;
	final float[][] xold;int nextXold;//old positions before the collisions
	private final Box[] boxes;int nextBx;//these are used to re-cycle the box and point objects
	private final Point[] points;int nextPt;
	private final Point[] ptXref; //point to the position in points that represents the point that is added because of the particle that is at position i in the index ArrayList ptrs 
	final private ArrayListInt NNList=new ArrayListInt(120);
	private ArrayListInt ptrs;//maps the offsets to global positions in the x,y,z vectors

	public static KDTree MMACreateEngine(double[] cx,double[] cy,double[] cz,int[] ints){
		ArrayListInt ptr=new ArrayListInt(ints);
		KDTree eng=new KDTree(ptr.length,cx,cy,cz);
		eng.rebuild(ptr);
		return eng;
	}
	public Point getPt(int off){return ptXref[off];}

	/**
	 * create a KDTree that will try to use the positions stored in the main system arrays. The KDTree created by this method is empty and will only be filled when a call to build() is called
	 * @param nn number of pk/cell you expect
	 * @param cx, cy,cz array's that have the position information stored
	 */
	KDTree(int nn,final double[] cx,final double[] cy,final double[] cz){
		x=cx;y=cy;z=cz;
		ptXref=new Point[conservativeUpperBound(nn)];
		nn*=2;//effectively because all particles may collide in one step this is practically an uppper bound on the number of points that will live in a cell.
		boxes=new Box[conservativeUpperBound(2*nn)];//upper bound on number of boxes that will be needed to store n points
		for(int i=0;i<boxes.length;i++)boxes[i]=new Box();
		points=new Point[conservativeUpperBound(3*boxes.length)];////each box has 3 points in it
		for(int i=0;i<points.length;i++)points[i]=new Point();
		xold=new float[points.length][3];

		Point.setStaticVariables(xold,x,y,z);
	}

	public KDTree(Gas gas, Parameters pars) {
		this(pars.pkPerCell,gas.Cx,gas.Cy,gas.Cz);
	}

	private void reset(){
		nextBx=nextPt=nextXold=0;//reset where we point
		final Point ptLL=popPt(),ptHH=popPt();//upper and lower bounds of the complete points
		setPt(ptLL,-BIG,-BIG,-BIG);
		setPt(ptHH,BIG,BIG,BIG);
		root=popBox();
		root.set(null,null,null,ptLL,ptHH,null);
	}
	final private static int conservativeUpperBound(double size){return (int)ceil(size+6*sqrt(size));}

	final private void setPt(final Point pt,final float x0,final float x1,final float x2){
		pt.set(-1, nextXold);
		nextXold++;
		pt.set(x0, x1, x2);
	}
	final private void setPt(final Point pt,final float x0,final float x1,final float x2,final int dim,final float cx){//this is hidious but It's all done in the name of efficiancy
		if(dim==0)     setPt(pt,cx,x1,x2);
		else if(dim==1)setPt(pt,x0,cx,x2);
		else		   setPt(pt,x0,x1,cx);
	}


	/** this method rebuilds the tree using the points pointed ty by ptrs array
	 * @param ptrs list of pointers to the points in the cell
	 */	public void rebuild(final ArrayListInt pkPtrs){
		 ptrs=pkPtrs;//update which cell pointer list we are using
		 reset();//unlink tree
		 //now we add all the particles:
		 for(int i=0;i<ptrs.length;i++)add(i);
	 }
	 /** add a point with a given offset
	  * @param off 	offset in global pkptrs
	  */
	 final private void add(int off){
		 final Point pt=popPt();
		 pt.set(ptrs.get(off), -1);//new point that refers to the particel's position
		 ptXref[off]=pt;//add Xref
		 //now we find/create a box that should contain this point:
		 Box parentBox=root;
		 int dim=0;//which dimension we are searching in
		 while(true){
			 if(parentBox.empty())break;
			 if(!parentBox.hasChildren())die(541,"bug: I'm supposed to have childern!");//TODO delme debugging
			 if(pt.x(dim)<parentBox.ctr.x(dim))parentBox=parentBox.dau1;//lower box
			 else parentBox=parentBox.dau2;//higher box
			 dim=nextDim(dim);//
		 }
		 addPtToBox(parentBox,pt,dim);
	 }



	 private void addPtToBox(final Box box,final Point pt,final int dim) {
		 if(!box.empty())die(999,"bug man I hat this!");//TODO: DEBUG DELME
		 box.ctr=pt;
		 //split me and create my daughters;
		 final Point dau1hi=popPt(),dau2lo=popPt();//now we need to set these pks
		 final float cx=box.ctr.x(dim);//the dimention we are splitting over
		 final Point lo=box.lo,hi=box.hi;
		 setPt(dau1hi,hi.x(0),hi.x(1),hi.x(2),dim,cx);
		 setPt(dau2lo,lo.x(0),lo.x(1),lo.x(2),dim,cx);
		 if(box.dau1!=null||box.dau2!=null)die(1);//TODO:debugging delme eventually. 
		 box.dau1=popBox();box.dau2=popBox();
		 box.dau1.set(box, null, null, lo, dau1hi, null);
		 box.dau2.set(box, null, null, dau2lo,hi,null);
	 }

	 /** return a new box but with no guarantees on the state of the box*/
	 private Box popBox(){
		 //we only create objects when they are needed...
		 final Box res=boxes[nextBx];
		 nextBx++;
		 return res;
	 }

	 /** return a new Point but with no guarantees on the state of the Point*/
	 final private Point popPt(){
		 //we only create objects when they are needed...
		 final Point res=points[nextPt];
		 nextPt++;
		 return res;
	 }

	 /**make a point's position be stored in the local cache insted of the global x,y,z take care of finding a pointer to it*/ 
	 public void cache(final int off) {
		 final Point pt=ptXref[off];//point in question
		 final float xf=pt.x(0),yf=pt.x(1),zf=pt.x(2);//yes, I'm not using Point's get method "premature optimization....."
		 setPt(pt,xf,yf,zf);
	 }
	 /** return a list of NN based on this KDTree. Only return the final position of points. 
	  * 
	  * @param offJ offsent in pkPointers of the point
	  * @param eps distance 
	  * @return
	  */
	 //don't want to re-initialize this every call to NN
	 static Box[] task=new Box[50];//should be enough for ~O(2^50) points.... I think

	 public ArrayListInt NN(final int offJ,final double eps) {
		 final double eps2=eps*eps;
		 if(eps<0)die(123211,"dude, comeon give me a positive eps damn in!");//TODO:Delme debug only 
		 NNList.reset();
		 int dim=0;

		 final Point pt=ptXref[offJ]; 
		 //this code is slightly modified from the NR code. So yes I don't really know how it works :) (see page 1110 SS 21.2)
		 Box enclosingBox=root;
		 //this block finds the smallest box that contains a spher of radious eps 
		 while(enclosingBox.hasChildren()){
			 Box oldBox=enclosingBox;
			 Box d1=enclosingBox.dau1;
			 Box d2=enclosingBox.dau2;
			 if(pt.x(dim)+eps<=d1.hi.x(dim))enclosingBox=d1;//if higher point is in the box the lower point is also in the box
			 else if (pt.x(dim)-eps>=d2.lo.x(dim)) enclosingBox=d2;//same idea here
			 dim=nextDim(dim);
			 if(enclosingBox==oldBox)break;//no change in box means that current box is smallest box
		 }

		 task[1]=enclosingBox;
		 int ntask=1;//number of pending jobs
		 while(ntask!=0){
			 Box k=task[ntask--];
			 if(k.hasChildren()&&k.ctr==null)die(411,"this is .... impossible!");//TODO: delme debug sanity check ...
			 if(Box.dist2(k,pt)>eps2)continue;//I'm too far away  
			 if(k.ctr==null)continue;//k is childless and pointless... like many people I know
			 if(!k.ctr.cached()&&Point.dist2(k.ctr,pt)<eps2)
				 NNList.add(k.ctr.globalRef);
			 task[++ntask]=k.dau1;
			 task[++ntask]=k.dau2;

		 }
		 // TODO Auto-generated method stub
		 return NNList;
	 }
	 /** update the position of particle pk to it's new location.
	  * this method will automatically clone the point in question
	  * @param off	offset of point in question in pkPtrs
	  */
	 public void update(final int off) {
		 //find the Point object that needs to be updated
		 final Point pt=ptXref[off];//Point object that represents the particle//TODO: DELME NOT NEEDED OR USED.
		 //sanity check::confirm that it is not cached://TODO DELME DEBUG
		 if(!pt.cached())die(321,"I AM supposed to be cached");//TODO:DELME DEBUG
		 //add the point again to the tree
		 add(off);
		 //now we have 2 points that originally come from the same particle (obviously at different times) 
	 }

	 final static private int nextDim(int dim) {return ++dim%3;}// I just want this to be inlined :)

	 public String toString(){
		 return "Graphics3D[{"+root.toString()+"},Axes->True]";
	 }

}



/** this is my point object. this thing should never explicitly hold some coordinates
 * @author husain
 */
class Point{
	private static double[] x,y,z;//this variable has a reference to the global postion vectors in the main gas
	private static float[][] xcache;//position of points that are not stored in the main gas
	public int globalRef=-1;//TODO: delete the public modifier 
	public int cacheRef=-1;//TODO: delete the public modifier 

	final public float x(int dim) {
		if(globalRef>-1){
			if(dim==0)return (float)x[globalRef];
			if(dim==1)return (float)y[globalRef];
			return (float)z[globalRef];
		}
		return xcache[cacheRef][dim];
	}

	public void set(float x0, float x1, float x2) {
		if((globalRef>-1)||(cacheRef<0))die(999);//TODO: delme debug make sure we don't touch a point that was not cloned
		xcache[cacheRef][0]=x0;
		xcache[cacheRef][1]=x1;
		xcache[cacheRef][2]=x2;
	}

	/** change the current point to be identical to the argument pt
	 * @param pt	argument pt
	 */
	public void set(final int gr,final int cr) {this.globalRef=gr;this.cacheRef=cr;//ta da!
	}

	public boolean cached(){return (cacheRef>-1);}

	public static void setStaticVariables(final float[][] cached2D,final double[]xg,final double[] yg,final double[]zg) {
		if(xcache!=null)die(4,"hmm were we called twise? ");//TODO: delme debug
		x=xg;y=yg;z=zg;
		xcache=cached2D;
	}

	/** reset what this Point represents. make both pointers -1
	 */
	public void reset(){globalRef=cacheRef=-1;}
	final static float dist2(final Point p,final Point q){
		return sq( q.x(0)-p.x(0))+ sq( q.x(1)-p.x(1))+ sq( q.x(2)-p.x(2));
	}
	public String toString(){
		final float LG=6*450;
		if(max(globalRef,cacheRef)<0)return "{gR&cR<0}";//TODO:delme debug
		float a=this.x(0),b=this.x(1),c=this.x(2);
		if(a>LG)a=LG;if(a<-LG)a=-LG;
		if(b>LG)b=LG;if(b<-LG)b=-LG;
		if(c>LG)c=LG;if(c<-LG)c=-LG;
		return "{"+a+","+b+","+c+"}";
	}
}

/** this is the main class for the boxes that populate the KD-Tree
 */
class Box{
	Point lo,hi,ctr;//these are my edgies;; ctr is the center point (ie. payload if it exists)
	Box mom,dau1,dau2;//parant node & it's children:: daughter 2 is older than #1 along dimention splitDim


	public void set(Box mum,Box d1,Box d2,Point l,Point h,Point c){mom=mum;dau1=d1;dau2=d2;lo=l;hi=h;ctr=c;}

	public void set(Box mom2, Point lo2, Point hi2) {mom=mom2;lo=lo2;hi=hi2;//create an childless box
	ctr=null;dau1=dau2=null;
	}


	boolean empty(){return ctr==null;}//we have an empty box
	boolean hasChildren(){return dau1!=null;}
	final static float dist2(final Box b,final Point p){
		float dd=0;
		if(p.x(0)<b.lo.x(0))dd+=sq(p.x(0)-b.lo.x(0));if(p.x(0)>b.hi.x(0))dd+=sq(p.x(0)-b.hi.x(0));
		if(p.x(1)<b.lo.x(1))dd+=sq(p.x(1)-b.lo.x(1));if(p.x(1)>b.hi.x(1))dd+=sq(p.x(1)-b.hi.x(1));
		if(p.x(2)<b.lo.x(2))dd+=sq(p.x(2)-b.lo.x(2));if(p.x(2)>b.hi.x(2))dd+=sq(p.x(2)-b.hi.x(2));
		return dd;
	}

	public String toString(){
		String res="";
		res="{Opacity[0.0],EdgeForm[Hue["+nextDouble()+"]],Cuboid["+lo.toString()+","+hi.toString()+"]}";
		if(!this.empty())res+=", {PointSize[.0175], Point["+ctr.toString()+"]}";
		if(this.hasChildren()){
			res+=","+this.dau1.toString()+","+this.dau2.toString();
		}
		return res;
	}

}

