/** this code is a modified version of the NN routine in the Numerical Recipes Book
 *  this code is not easy to understand because (as much as I adore NR) they write their code in a non-understandable way. 
 *  I have tried to modify things but only slightly. I did not do a complete re-write because:
 *  1. I'm lazy and
 *  2. it actually may be more efficant 
 *  I will try to remove the Point class from this implementation and directly manage gas.cx,cy,cz.
 */

package edu.mit.husain.DSMCVR3D;

import java.io.Serializable;



public class SlowKDtree implements Serializable{


	/**
	 * 
	 */
	private static final long serialVersionUID = 3L;
	final static int DUMMY=100;
	static final double BIG=1.0e99;//edge of boxes
	static int nboxes, npts; //number of boxes and number of points later should be changed....
	public static SlowPoint[] ptss;//coordinates of all points
	static Boxnode[] boxes=new Boxnode[DUMMY];

	static int[] ptindx=new int[DUMMY],rptindx=new int[DUMMY];
	static double[] coords=new double[DUMMY];


	static int cp_tdimnpts;//stupid variables to fake pointer artheithmetic -- OO indeed...
	static int hp_ptlo;//same thing.....

	/** Default Constructor
	 * 
	 */
	SlowKDtree(){
		//if needed
	}

	/** This is the constructor of the NN Engine. this is all very slow ...  
	 * @param x Cx
	 * @param y Cy
	 * @param z Cz
	 * @param ptrs Integer ArrayList of the pointers to all the particles in a cell. 
	 * @author husain + NR
	 */

	SlowKDtree(final double[] x,final double[] y,final double[] z,ArrayListInt cellPtrs){
		//I don't know if this is that useful. since this is so slow I might as well 
		System.gc();//clean garbage ...

		ptss=null;
		System.gc();//this is mess, i might as well delete things.
		ptss=new SlowPoint[cellPtrs.length];
		for(int i=0;i<cellPtrs.length;i++){
			ptss[i]=new SlowPoint(x[cellPtrs.get(i)],y[cellPtrs.get(i)],z[cellPtrs.get(i)]);
		}
		npts=cellPtrs.length;
		if(ptindx.length!=npts){ptindx=new int[npts];rptindx=new int[npts];}

		int ntmp,m,k,kk,j,nowtask,jbox,np,tmom,tdim,ptlo,pthi;
		//
		int[] taskmom=new int[50],taskdim=new int[50];
		for(k=0;k<npts;k++)ptindx[k]=k;//initialize index of monts


		//calculate number of boxes:
		m=1;
		for(ntmp=npts;ntmp!=0;ntmp>>=1){//god i hate this 
			m<<=1;
		}
		nboxes=2*npts-(m>>1);
		if(m<nboxes)nboxes=m;
		nboxes--;

		if(boxes.length!=nboxes){ boxes=new Boxnode[nboxes];}//
		//copy points into coords.... if u have an arryay of points why copy it again ??!!! :(
		if(coords.length!=3*npts){ coords=new double[3*npts];}
		for(j=0,kk=0;j<3;j++,kk+=npts){
			for(k=0;k<npts;k++)coords[kk+k]=ptss[k].x[j];
		}
		//do the boxes now ...
		SlowPoint lo=new SlowPoint(-BIG,-BIG,-BIG),hi=new SlowPoint(BIG,BIG,BIG);
		boxes[0]=new Boxnode(lo,hi,0,0,0,0,npts-1);//root box
		jbox=0;
		taskmom[1]=0;

		taskdim[1]=0;
		nowtask=1;
		while(nowtask!=0){

			tmom=taskmom[nowtask];
			tdim=taskdim[nowtask--];
			ptlo=boxes[tmom].ptlo;
			pthi=boxes[tmom].pthi;
			hp_ptlo=ptlo;//hp= &ptindx[ptlo];

			cp_tdimnpts=tdim*npts;//cp= &coords[tdim*npts];
			np=pthi-ptlo+1;
			kk=(np-1)/2;
			selecti(kk,np);//selecti(kk,hp,np,cp);
			hi.assign(boxes[tmom].hi);
			lo.assign(boxes[tmom].lo);
			hi.x[tdim]=lo.x[tdim]=coords[tdim*npts+indx(kk)];//coords[tdim*npts+hp[kk]];
			boxes[++jbox]=new Boxnode(boxes[tmom].lo,hi,tmom,0,0,ptlo,ptlo+kk);
			boxes[++jbox]=new Boxnode(lo,boxes[tmom].hi,tmom,0,0,ptlo+kk+1,pthi);
			boxes[tmom].dau1=jbox-1;
			boxes[tmom].dau2=jbox;

			if(kk>1){
				taskmom[++nowtask]=jbox-1;
				taskdim[nowtask]=(tdim+1)%3;
			}
			if(np-kk>3){

				taskmom[++nowtask]=jbox;
				taskdim[nowtask]=(tdim+1)%3;
			}
		}
		for(j=0;j<npts;j++) rptindx[ptindx[j]]=j;

	}
	public 	SlowKDtree(final double[] x,final double[] y,final double[] z,int[] cellPtrs){
		this(x,y,z,new ArrayListInt(cellPtrs));
	}

	public int[] MMAbenchmark(final double[]cx,final double[] cy,final double[] cz,int[] ptrsInt,int[] searchPk,double eps){
		ArrayListInt ptrs=new ArrayListInt(ptrsInt);
		ArrayListInt res=new ArrayListInt(10);
		SlowPoint pt=new SlowPoint(0,0,0);
		for(int i:searchPk){
			pt.assignVsAtPos(cx, cy, cz, i);
			locatenear(pt,eps,res,ptrs);
		}
		return res.trimmedArray();
	}

	public int locatenear(final SlowPoint pt,final double r,ArrayListInt NNList,ArrayListInt ptrs){
		NNList.reset();
		
		int k,i,nb,nbold,ntask,jdim,d1,d2;
		int[] task=new int[50];
		nb=jdim=0;
		if(r<0.0)System.err.println("r should be >0");
		while(boxes[nb].dau1!=0){
			nbold=nb;
			d1=boxes[nb].dau1;
			d2=boxes[nb].dau2;
			if(pt.x[jdim]+r<=boxes[d1].hi.x[jdim])nb=d1;
			else if (pt.x[jdim]-r>=boxes[d2].lo.x[jdim]) nb=d2;
			jdim=++jdim%3;//come-on this just ugly ....
			if(nb==nbold)break;
		}
		task[1]=nb;

		ntask=1;
		while(ntask!=0){
			k=task[ntask--];

			if(dist2(boxes[k],pt)>r*r)continue; 
			if(boxes[k].dau1!=0){

				task[++ntask]=boxes[k].dau1;
				task[++ntask]=boxes[k].dau2;
			}else{
				for(i=boxes[k].ptlo;i<=boxes[k].pthi;i++){
					if((dist2(ptss[ptindx[i]],pt)<=r*r))
						NNList.add(ptrs.get(ptindx[i]));
				}
			}
		}
		return NNList.length;
	}

	

	
	static double dist2(final SlowBox b,final SlowPoint p){
		double dd=0;
		for(int i=0;i<3;i++){
			if(p.x[i]<b.lo.x[i])dd+=SQR(p.x[i]-b.lo.x[i]);
			if(p.x[i]>b.hi.x[i])dd+=SQR(p.x[i]-b.hi.x[i]);
		}
		return dd;
	}

	static double dist2(final SlowPoint p,final SlowPoint q){
		double dd=0.0;
		for(int j=0;j<3;j++)dd+= SQR( q.x[j]-p.x[j]);
		return dd;
	}



	final static double SQR(double x){return x*x;}

	int selecti(final int k,int n){//p1104
		//arr is not modified......
		int i,ia,ir,j,l,mid;
		double a;


		l=0;
		ir=n-1;//not clear to me if n is just the length of the indx variables
		for(;;){
			if(ir<=l+1){

				if(ir==l+1&&arr(indx(ir))<arr(indx(l))){SWAPindx(l,ir);}
				return indx(k);

			}else{
				mid=(l+ir)>>1;//fancy div by 2
				SWAPindx(mid,l+1);
				if(arr(indx(l))>arr(indx(ir)))SWAPindx(l,ir);
				if(arr(indx(l+1))>arr(indx(ir)))SWAPindx(l+1,ir);
				if(arr(indx(l))>arr(indx(l+1)))SWAPindx(l,l+1);
				i=l+1;

				j=ir;
				ia=indx(l+1);
				a=arr(ia);

				for(;;){
					do i++; while(arr(indx(i))<a);
					do j--; while(arr(indx(j))>a);
					if(j<i)break;
					SWAPindx(i,j);

				}
				indxeq(l+1,indx(j));
				indxeq(j,ia);//indx[j]=ia;

				if(j>=k)ir=j-1;
				if(j<=k)l=i;
			}

		}
	}
	int indx(int i){return ptindx[i+hp_ptlo];}//i hate u c++
	void SWAPindx(int i,int j){SWAP(ptindx,i+hp_ptlo,j+hp_ptlo);}
	void indxeq(int i,int x){ptindx[i+hp_ptlo]=x;}
	double arr(int i){return coords[i+cp_tdimnpts];}



	static void SWAP(int[] a,int i,int j){//swap i and j
		final int t=a[i];
		a[i]=a[j];a[j]=t;
	}

	static void assign(int[] a,final int x){for(int i=0;i<a.length;i++)a[i]=x;}}

class SlowBox{
	SlowPoint lo,hi;
	SlowBox(){}
	SlowBox(SlowPoint mylow,SlowPoint myhi){
		lo=new SlowPoint(mylow);//it's hard to figure out what C++ constructors do in the NR book...
		hi=new SlowPoint(myhi);
	}
}

class Boxnode extends SlowBox{
	int mom,dau1,dau2,ptlo,pthi;
	Boxnode(SlowPoint mylo,SlowPoint myhi,int mymom,int myd1,int myd2,int myptlo,int mypthi){
		super(mylo,myhi);
		mom=mymom;dau1=myd1;dau2=myd2;ptlo=myptlo;pthi=mypthi;
	}
}

class SlowPoint{

	final double[] x;
	SlowPoint(double x1,double y,double z){x=new double[3];
	x[0]=x1;
	x[1]=y;
	x[2]=z;
	}
	SlowPoint(SlowPoint pt){this(pt.x[0],pt.x[1],pt.x[2]);}
	void assign(SlowPoint pt){x[0]=pt.x[0];x[1]=pt.x[1];x[2]=pt.x[2];}
	public void assignVs(double xx,double yy,double zz){x[0]=xx;x[1]=yy;x[2]=zz;}
	public String toString(){return "{"+x[0]+","+x[1]+","+x[2]+"}";}
	public void assignVsAtPos(double[] cx, double[] cy, double[] cz, int i) {
		cx[i]=x[0];cy[i]=x[1];cz[i]=x[2];
	}
	
	public double[] toArray(){return x;}//used for debuging

}
