// XYZbins.h
// Robert Baertsch 
// 28 Feb 2002
//
// Modified by Kevin Karplus
//	8 April 2002
//	removed lots of extraneous stuff (some of which belongs
//		in specific programs, not in library)
//	removed timing data

#ifndef XYZBINS_H
#define XYZBINS_H

#include "Pointlist.h" 
#include "IntArray.h"
#include <assert.h>

class XYZbins
{
  private:
    
    const Pointlist& Points;	// what points are in the bins
    // Points is NOT owned by XYZbins, so whoever creates
    // XYZbins is responsible for making sure that the Pointlist
    // doesn't change while the XYZbins are around.
    
    // Each bin is a list of the subscripts (in Points) of
    // of the points that fall within the cube defined by the bin.
    //
    // get_subscripts() will compute the 3 subscripts of the bin
    // for a particular point.
    //
    // subscript() will conver the 3 subscripts to a single subcript
    // for the Bins array.
    IntArray **Bins;
    
    double ScaleX ; // multiplier to calc which bin atom is in (x coord)
    double ScaleY ; // multiplier to calc which bin atom is in (y coord)
    double ScaleZ ; // multiplier to calc which bin atom is in (z coord)
    
    int NumBinsX;  // number of bins along x axis
    int NumBinsY;  // number of bins along y axis
    int NumBinsZ;  // number of bins along z axis
    
    XYZpoint minp ;   // smallest x,y,z coordinates of protein  
    XYZpoint maxp ;   // largest x,y,z coordinates of protein  


    inline int subscript(int i, int j, int k)
    {    return ((i*NumBinsY)+j)*NumBinsZ+k;
    }
  
    // compute a subscript for one dimension.
    inline int get_one_subscript(double value,
    	double min_value,
	double max_value,
	double scale,
	int num_bins)
    {   if (value<= min_value) return 0;
    	if (value>= max_value) return num_bins-1;
	int sub = static_cast<int>(floor(scale*(value-min_value)));
	assert(sub >=0 && sub < num_bins);
	return sub;
    }
    
    inline void get_subscripts(int &i, int &j, int &k, 
    	const XYZpoint& A)
    {   if (NumBinsX==0) 
    	{   i=j=k=0; 	// empty bins
	    return;
	}
	i= get_one_subscript(A.X, minp.X, maxp.X, ScaleX, NumBinsX);
    	j= get_one_subscript(A.Y, minp.Y, maxp.Y, ScaleY, NumBinsY);
    	k= get_one_subscript(A.Z, minp.Z, maxp.Z, ScaleZ, NumBinsZ);
    }
  
  private:
    IntArray*& bin_element(int i, int j, int k)
    {	return Bins[subscript(i,j,k)];
    }

    IntArray*& bin_element(const XYZpoint& A)
    {	int i,j,k;
        get_subscripts(i,j,k,A);
        return Bins[subscript(i,j,k)];
    }

    // Report the subscript of the closest point to spot
    // counting only points in the same bin as spot.
    // Returns -1 if no point in same bin.
    //
    // dist2 is set to sq_dist(spot, Points[return value]), if point found.
    int approx_closest(double &dist2, const XYZpoint& spot);

  public:
    // This constructor creates the bins and fills them with 
    // all the points in the subarray of Atoms.
    //
    // Points with any coordinates >= MaxCoord are omitted.
    //	(MaxCoord defined in XYZbins.cc as     const double MaxCoord= 9.e6;
    //
    //
    // As in Pointlist, fix_start_length() is used to interpret
    // negative values of start and length:
    // start<0 is interpreted as
    //    Atoms.num_points()+start (so the last element of the list is -1).
    // Length<0 is interpreted in a similar way, so that
    //    (start,-1) means from start to the end of the list,
    //    (start,-2) omits the last element, and so on.
    //
    // If atom_present is provided, then only those points whose
    //	atom_present[p] is true will be included.
    // The subscripting is identical in Atoms and atom_present
    
    XYZbins(const Pointlist& Atoms, 
    	double BinWidth=8.0,
	const bool *atom_present=NULL,
	int start=0, int length=-1
    	);


    // This constructor fills the bins only with selected points---
    // those whose subscripts are in
    // selected_atoms[0]..selected_atoms[num_selected-1] 
    XYZbins(const Pointlist& Atoms, 
	const int *selected_atoms, 
	int num_selected,
	double BinWidth=8.0
    	);


    ~XYZbins(void)
    {   for (int i=NumBinsX*NumBinsY*NumBinsZ-1; i>=0; i--)
	{    if (Bins[i]!= 0)
		delete Bins[i];
	}
	delete [] Bins;
    }

    // There are 3 closely related routines for counting the number
    // of points within a fixed radius of a given spot.
    // All are called num_points_near().
    // One takes a single query point, one a Pointlist of queries,
    // and one (which shouldn't be used in new code) 2 Pointlists of queries.
    
    // Counts points in the bins that are within spot_radius
    //    Angstroms of point spot 
    int num_points_near(const XYZpoint& spot,  double spot_radius);

    // Query the bins data structure repeatedly
    // storing the result of num_points_near for each spot in pre-allocated
    // array burial.
    // As in Pointlist, fix_start_length() is used to interpret
    // negative values of start and length:
    // The burial[0] value is for spot (*spots)[start].
    // 
    // Call is a no-op if spots==NULL.
    void num_points_near(unsigned short* burial,
        const Pointlist *spots,
        double spot_radius,
        int start=0, int length=-1 
	);

    // Query repeatedly, interleaving the queries for spots1 and spots2.
    // Preallocated burial array should be TWICE the length.
    // NOTE: use of Query2ListSpots is deprecated---it is
    // provided only for compatibility with old Kestrel-based code.
    //
    // If spots1==NULL or spots2==NULL, this code is equivalent
    // to QueryListSpots---no interleaving is done.
    void num_points_near(unsigned short* burial,
        const Pointlist *spots1,
        double spot_radius1,
        const Pointlist *spots2,
        double spot_radius2,
        int start=0, int length=-1 
	);


    // Similar to the num_points_near() members,
    // there are two apply_to_points_near() members, 
    // which apply a function to each pair of close points found.
     
    typedef void ClosePairFcn(int index_in_bins, const XYZpoint& bin_point,
    			int in_query, const XYZpoint& query_point,
			double dist2);
    
    // apply the function f(index_in_bins, bin_point, spot_index, spot, dist^2)
    // to each bin_point within spot_radius of spot.
    void apply_to_points_near(ClosePairFcn *f,
    	const XYZpoint& spot,  double spot_radius,
	int spot_index=1);


    // apply the function f(index_in_bins, bin_point, i, spots[i], dist^2)
    // to each bin_point within spot_radius of spots[i].
    void apply_to_points_near(ClosePairFcn *f,
        const Pointlist *spots,
        double spot_radius,
        int start=0, int length=-1 
	);



    // Report the subscript of the closest point to spot
    // counting only points with squared distance <= max_dist2.
    // (Unless max_dist2< 0, then put no constraints on search.)
    //
    // Returns -1 if no point meeting constraints can be found.
    //
    // dist2 is set to sq_dist(spot, Points[return value]), if point found.
    int closest(double &dist2, const XYZpoint& spot, double max_dist2=-1.);
    
    // Return the subscript of the closest point to any points in 
    // subarray of spots (specified by start and length).
    // counting only pairs with squared distance <= max_dist2.
    // (Unless max_dist2< 0, then put no constraints on search.)
    // Returns -1 if no pair meeting constraints can be found.
    // spot_subscript is set to the index in spots of the point closest
    // to the returned point.
    // (Won't report spots for which spot_present[s]=false, if
    //    spot_present provided.)
    int closest(double &dist2,
		int & spot_subscript,
		const Pointlist* spots,
		double max_dist2=-1.,
		const bool *spot_present=NULL,
		int start=0, int length=-1);

};


// CHANGE LOG:
// 18 April 2002 Kevin Karplus
//	Added new constructor for selected subset of points.
// 27 Feb 2004 Kevin Karplus
//	Updated old-style casts
// Tue Jul 26 12:42:05 PDT 2005 Kevin Karplus
//	Added atom_present argument to basic constructor.
// Tue Jul 26 13:25:42 PDT 2005 Kevin Karplus
//	Added spot_present to closest 
#endif