// test_xyzbins.cc
//	Kevin Karplus
//	8 April 2002

// This is a crude test fixture for testing the XYZbins class.
//
// Two Pointlists are created:
//	one is used to set up XYZbins, the other to query the bins.
// The "correct" answers are determined by brute-force comparison of
// all points.

#include <stdlib.h>	// for drand48
#include "XYZbins.h"

const double XRange = 25. ;	// range of coordinates



XYZpoint gen_point(void)
{
  XYZpoint temp;
	
	temp.X = XRange*(drand48()-0.5);
	temp.Y = XRange*(drand48()-0.5);
	temp.Z = XRange*(drand48()-0.5);
	
  return temp;
}

static int Count_for_func;
void Count_reported_pairs(int index_in_bins, const XYZpoint& bin_point,
	int in_query, const XYZpoint& query_point,
	double dist2)
{
    Count_for_func++;
}

main()
{
    // atoms in dense part of proteins are about 40/cubic nanometer
    const int num_atoms = (int)ceil(XRange*XRange*XRange*0.040);
    const int num_queries = num_atoms/5;
    const int num_tries = 10;	// number of repetitions of test
    
    Pointlist atoms(num_atoms), queries(num_queries);
    XYZpoint A;
    for (int outer=0; outer<num_tries; outer++)
    {   // generate random lists
	atoms.clear();
	for(int i=0;i<num_atoms;i++)
	{   A = gen_point();
	    atoms.append_point(A);
	}
	queries.clear();
	for(int i=0;i<num_queries;i++)
	{   A = gen_point();
	    queries.append_point(A);
	}
	
	XYZbins bins(atoms);
	
	const double radius = 8.0;
	const double rad2 = radius*radius;
	
	cerr << "Checking burial\n";
	short unsigned int *burial = new
		short unsigned int[num_queries];
	bins.num_points_near(burial, &queries, radius);
	
	for(int i=0;i<num_queries;i++)
	{   int one_spot_burial = bins.num_points_near(queries[i], radius);
	    Count_for_func=0;
	    bins.apply_to_points_near(&Count_reported_pairs, queries[i], radius);
	    int real_burial = 0;
	    for(int a=0;a<num_atoms;a++)
	    {   if (sq_dist(queries[i], atoms[a]) <= rad2)
	    	    real_burial++;
	    }
	    if (burial[i] != real_burial  ||
	    	one_spot_burial != real_burial ||
		Count_for_func != real_burial)
	    {	cerr << "Error: Burial mismatch burial[" << i
	    	    << "]=" << burial[i]
		    << " one_spot_burial=" << one_spot_burial
		    << " function_applied=" << Count_for_func
	    		<< " real_burial=" << real_burial << "\n"
			<< "Spot "<< i<< " at " << queries[i] << "\n";
		
		for(int a=0;a<num_atoms;a++)
	    	{   if (sq_dist(queries[i], atoms[a]) <= rad2)
		    {    cerr << "Atom " << a << " at " << atoms[a]
		    		<< " sqdist="
				<< sq_dist(queries[i], atoms[a])
				<< "\n";
		    }
		}
		assert(burial[i]==real_burial);	// to get into debugger
	    }
	}
        delete [] burial;
	
	cerr << "Checking closest\n";
	for(int i=0;i<num_queries;i++)
	{   double dist2;
	    int closest_atom = bins.closest(dist2, queries[i]);
	    int real_closest = -1;
	    double max_dist2 = 1e99;
	    
	    for(int a=0;a<num_atoms;a++)
	    {   double dist2a = sq_dist(queries[i], atoms[a]);
	        if (dist2a <= max_dist2)
		{   real_closest=a;
		    max_dist2 = dist2a;
//		    cerr << "a=" << a << " dist2a=" << dist2a
//		    	<< " max_dist2=" << max_dist2
//			<< "\n";y
		}
	    }

	    assert(closest_atom==real_closest);
	    assert(max_dist2== dist2);

// WARNING: g++ version 2.91.66 and version 2.96
// both mess up the compilation of the test program with -O
// They seem to lose max_dist2 before getting to the assertion, 
// and set it to zero.

	}

    }
	

}