#include "Quaternion.h"
#include "XYZpoint.h"
#include "Transform.h"
#include "Pointlist.h"
#include <time.h>
#include <stdlib.h>

XYZpoint gen_point(void);

int main(void)
{
    XYZpoint A;
    long int seed, nump, num_tries;

    cout <<"\nGive me the number of points to transform ...";
    cin >> nump;
    cout<< "\nGive me a random seed ...";
    cin >> seed; 
    srand48(seed);
    
    cout<< "\nGive me the number of tries to average...";
    cin >> num_tries; 

    int superpose_iter = (5000+2*nump)/nump;
    int trans_iter=10*superpose_iter;
    int dist_iter = trans_iter;
    
    double  superpose_time =0;
    double  coplanar_time=0;
    double Transform_time  =0;
    double Dist2_time =0;
    Pointlist X(nump), Y(nump),  Z(nump);

    Transform T;
    for (int outer=0; outer<num_tries; outer++)
    {   
	X.clear();	Y.clear();
	for(int i=0;i<nump;i++)
	{
		A = gen_point();
		X.append_point(A);
		A = gen_point();
		Y.append_point(A);
	}
	
	double start_cputime = clock();	
	for (int i=1; i<superpose_iter; i++)
	{   T.optimal_transform(X,Y);
	}
	double new_cputime = clock();
	superpose_time += new_cputime - start_cputime;

	double start_trans_cputime = clock();
	for (int i=0; i<trans_iter; i++)
		T.trans(X,Z);
	new_cputime = clock();
	Transform_time  += new_cputime - start_trans_cputime;
	
	double start_square_time = clock();
	for (int i=0; i<dist_iter; i++)
		 square_dist(Z,Y);
	new_cputime = clock();
	Dist2_time += new_cputime - start_square_time;

	if (nump==3)
	{   double cputime = clock();	
	    for (int i=1; i<superpose_iter; i++)
	    {   T.coplanar_trans(X,Y);
	    }
	    new_cputime = clock();
	    coplanar_time += new_cputime - cputime;
	}
    }
    cout << "\n\n"
    	<< "For " << nump << " points, "
    	 << " optimal superpose (with eigenvectors)= " <<superpose_time/(superpose_iter*num_tries)
	 << " (" << superpose_time/(superpose_iter*num_tries*nump)
	 <<	" per point)\n"
	 << " transform/point= " << Transform_time/(num_tries*nump*trans_iter)
	 << " distance^2 time= " << Dist2_time/(num_tries*dist_iter)
	 << " microseconds\n";
    cout << "Combined superpose, transform, dist^2= "
    	<< (superpose_time/superpose_iter + Transform_time/trans_iter
		+Dist2_time/dist_iter) / num_tries
	<< " microseconds\n\n";
    if (nump==3)
    {    cout << "Coplanar transform time= " <<coplanar_time/(superpose_iter*num_tries)
	 << " microseconds, (" << coplanar_time/(superpose_iter*num_tries*nump)
	 <<	" per point)\n";
    }
    return 0;
}

XYZpoint gen_point(void)
{
  XYZpoint temp;
	
	temp.X = 100*drand48();
	temp.Y = 100*drand48();
	temp.Z = 100*drand48();
	
  return temp;
}

// 8 June 2004 Kevin Karplus
//	Made the "main" procedure have legal int return.
//	Increased information output (dist^2 time)