// test_array.cc -- test driver for array.hh classes
// $Id$

#include <cassert>
#include <math.h>
#include <cstdio>
using std::printf;

#include "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"

using jtutil::how_many_in_range;
using jtutil::fuzzy;
using jtutil::array1d;
using jtutil::array2d;
using jtutil::array3d;
using jtutil::array4d;

//
// prototypes for non-class functions defined in this file
//
void test_array1d_double(bool print_flag);
void test_array2d_float(bool print_flag);
void test_array3d_double(bool print_flag);
void test_array4d_double(bool print_flag);

//******************************************************************************

//
// Usage:
//	test_array  [ print ]
// If an argument is supplied (it doesn't matter what it is), then
// the array indices will be printed as they're tested.
//
int main(int argc, const char* const argv[])
{
const bool print_flag = (argc > 1);
test_array1d_double(print_flag);
test_array2d_float(print_flag);
test_array3d_double(print_flag);
test_array4d_double(print_flag);
}

//******************************************************************************

//
// This function tests the  array1d<double>  class.
//
void test_array1d_double(bool print_flag)
{
const int min_i = 11;		const int max_i = 13;
const double delta = 1.0/7.0;
const double change_value = 144.025036146;

printf("\ntesting array1d<double>...\n");

printf("constructing...\n");
array1d<double> A(min_i, max_i);
assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
assert(A.N_i() == how_many_in_range(min_i,max_i));

printf("assigning values...\n");
#define VALUE_1D(ix)	\
		double(1000.0 + ix + delta)
	for (int i = min_i ; i <= max_i ; ++i)
	{
	if (print_flag)
	   then printf("   assigning %d\n", i);
	A(i) = VALUE_1D(i);
	}

printf("checking N-D vs 1-D addresses...\n");
int sub = 0;
const double* array_ptr = A.get_array();
	for (int i = min_i ; i <= max_i ; ++i)
	{
	if (print_flag)
	   then printf("   checking %d\n", i);
	assert( &A(i) == &array_ptr[sub++] );
	}
assert(sub == A.N_array());

printf("checking which things change when we assign an array element...\n");
double sumsq_of_diff = 0.0;
	for (int i = min_i ; i <= max_i ; ++i)
	{
	if (print_flag)
	   then printf("   changing %d\n", i);
	double save_a = A(i);
	A(i) = change_value;

		for (int ii = min_i ; ii <= max_i ; ++ii)
		{
		if (print_flag)
		   then printf("      checking %d\n", ii);

		double should_be = (ii == i) ? change_value : VALUE_1D(ii);
		assert( fuzzy<double>::EQ(A(ii), should_be) );

		double diff = A(ii) - should_be;
		sumsq_of_diff += diff*diff;
		}
	A(i) = save_a;
	}
double RMS_diff = sqrt(sumsq_of_diff / A.N_array());

printf("==> everything looks ok (RMS_diff=%.3g)\n", RMS_diff);
}

//******************************************************************************

//
// This function tests the  array2d<float>  class.
//
void test_array2d_float(bool print_flag)
{
const int min_i = 11;		const int max_i = 13;
const int min_j = 21;		const int max_j = 24;
const float delta = 1.0/7.0;
const float change_value = 144.025036146;

printf("\ntesting array2d<float>...\n");

printf("constructing...\n");
array2d<float> A(min_i,max_i, min_j,max_j);
assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
assert(A.N_i() == how_many_in_range(min_i,max_i));
assert(A.N_j() == how_many_in_range(min_j,max_j));

printf("assigning values...\n");
#define VALUE_2D(ix,jx)	\
		float(1000.0 + 10.0*jx + ix + delta)
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	if (print_flag)
	   then printf("   assigning %d %d\n", i, j);
	A(i,j) = VALUE_2D(i,j);
	}
	}

printf("checking N-D vs 1-D addresses...\n");
int sub = 0;
const float* array_ptr = A.get_array();
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	if (print_flag)
	   then printf("   checking %d %d\n", i, j);
	assert( &A(i,j) == &array_ptr[sub++] );
	}
	}
assert(sub == A.N_array());

printf("checking which things change when we assign an array element...\n");
double sumsq_of_diff = 0.0;
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	if (print_flag)
	   then printf("   changing %d %d\n", i, j);
	float save_a = A(i,j);
	A(i,j) = change_value;

		for (int ii = min_i ; ii <= max_i ; ++ii)
		{
		for (int jj = min_j ; jj <= max_j ; ++jj)
		{
		if (print_flag)
		   then printf("      checking %d %d\n", ii, jj);

		float should_be = ((ii == i) && (jj == j))
				  ? change_value : VALUE_2D(ii,jj);
		assert( fuzzy<float>::EQ(A(ii,jj), should_be) );

		double diff = A(ii,jj) - should_be;
		sumsq_of_diff += diff*diff;
		}
		}
	A(i,j) = save_a;
	}
	}
float RMS_diff = sqrt(sumsq_of_diff / A.N_array());

printf("==> everything looks ok (RMS_diff=%.3g)\n", double(RMS_diff));
}

//******************************************************************************

//
// This function tests the  array3d<double>  class.
//
void test_array3d_double(bool print_flag)
{
const int min_i = 11;		const int max_i = 13;
const int min_j = 21;		const int max_j = 24;
const int min_k = 31;		const int max_k = 35;
const double delta = 1.0/7.0;
const double change_value = 144.025036146;

printf("\ntesting array3d<double> with passed-in array...\n");

const int N_i = how_many_in_range(min_i,max_i);
const int N_j = how_many_in_range(min_j,max_j);
const int N_k = how_many_in_range(min_k,max_k);
const int N = N_i*N_j*N_k;
printf("allocating %d*%d*%d = %d underlying storage array...\n",
       N_i,N_j,N_k, N);
double *array = new double[N];

  // local scope for array
  {
printf("constructing array3d<double> object...\n");
array3d<double> A(min_i,max_i, min_j,max_j, min_k,max_k, array);
assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
assert(A.min_k() == min_k);	assert(A.max_k() == max_k);
assert(A.N_i() == how_many_in_range(min_i,max_i));
assert(A.N_j() == how_many_in_range(min_j,max_j));
assert(A.N_k() == how_many_in_range(min_k,max_k));

printf("assigning values...\n");
#define VALUE_3D(ix,jx,kx)	\
		double(1000.0 + 100.0*kx + 10.0*jx + ix + delta)
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	for (int k = min_k ; k <= max_k ; ++k)
	{
	if (print_flag)
	   then printf("   assigning %d %d %d\n", i, j, k);
	A(i,j,k) = VALUE_3D(i,j,k);
	}
	}
	}

printf("checking N-D vs 1-D addresses...\n");
int sub = 0;
const double* array_ptr = A.get_array();
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	for (int k = min_k ; k <= max_k ; ++k)
	{
	if (print_flag)
	   then printf("   checking %d %d %d\n", i, j, k);
	assert( &A(i,j,k) == &array_ptr[sub++] );
	}
	}
	}
assert(sub == A.N_array());

printf("checking which things change when we assign an array element...\n");
double sumsq_of_diff = 0.0;
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	for (int k = min_k ; k <= max_k ; ++k)
	{
	if (print_flag)
	   then printf("   changing %d %d %d\n", i, j, k);
	double save_a = A(i,j,k);
	A(i,j,k) = change_value;

		for (int ii = min_i ; ii <= max_i ; ++ii)
		{
		for (int jj = min_j ; jj <= max_j ; ++jj)
		{
		for (int kk = min_k ; kk <= max_k ; ++kk)
		{
		if (print_flag)
		   then printf("      checking %d %d %d\n", ii, jj, kk);

		double should_be = ((ii == i) && (jj == j) && (kk == k))
				   ? change_value : VALUE_3D(ii,jj,kk);
		assert( fuzzy<double>::EQ(A(ii,jj,kk), should_be) );

		double diff = A(ii,jj,kk) - should_be;
		sumsq_of_diff += diff*diff;
		}
		}
		}
	A(i,j,k) = save_a;
	}
	}
	}
double RMS_diff = sqrt(sumsq_of_diff / A.N_array());
printf("==> everything looks ok (RMS_diff=%.3g)\n", RMS_diff);

printf("exiting local scope ==> destroying array3d object...\n");
  }

printf("deleting underlying storage array...\n");
delete[] array;
}

//******************************************************************************

//
// This function tests the  array4d<double>  class.
//
void test_array4d_double(bool print_flag)
{
const int min_i = 11;		const int max_i = 13;
const int min_j = 21;		const int max_j = 24;
const int min_k = 31;		const int max_k = 35;
const int min_l = 41;		const int max_l = 46;
const double delta = 1.0/7.0;
const double change_value = 144.025036146;

printf("\ntesting array4d<double>...\n");

printf("constructing...\n");
array4d<double> A(min_i,max_i, min_j,max_j, min_k,max_k, min_l,max_l);
assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
assert(A.min_k() == min_k);	assert(A.max_k() == max_k);
assert(A.min_l() == min_l);	assert(A.max_l() == max_l);
assert(A.N_i() == how_many_in_range(min_i,max_i));
assert(A.N_j() == how_many_in_range(min_j,max_j));
assert(A.N_k() == how_many_in_range(min_k,max_k));
assert(A.N_l() == how_many_in_range(min_l,max_l));

printf("assigning values...\n");
#define VALUE_4D(ix,jx,kx,lx)	\
		double(10000.0 + 1000.0*kx + 100.0*jx + 10.0*ix + lx + delta)
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	for (int k = min_k ; k <= max_k ; ++k)
	{
	for (int l = min_l ; l <= max_l ; ++l)
	{
	if (print_flag)
	   then printf("   assigning %d %d %d %d\n", i, j, k, l);
	A(i,j,k,l) = VALUE_4D(i,j,k,l);
	}
	}
	}
	}

printf("checking N-D vs 1-D addresses...\n");
int sub = 0;
const double* array_ptr = A.get_array();
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	for (int k = min_k ; k <= max_k ; ++k)
	{
	for (int l = min_l ; l <= max_l ; ++l)
	{
	if (print_flag)
	   then printf("   checking %d %d %d %d\n", i, j, k, l);
	assert( &A(i,j,k,l) == &array_ptr[sub++] );
	}
	}
	}
	}
assert(sub == A.N_array());

printf("checking which things change when we assign an array element...\n");
double sumsq_of_diff = 0.0;
	for (int i = min_i ; i <= max_i ; ++i)
	{
	for (int j = min_j ; j <= max_j ; ++j)
	{
	for (int k = min_k ; k <= max_k ; ++k)
	{
	for (int l = min_l ; l <= max_l ; ++l)
	{
	if (print_flag)
	   then printf("   changing %d %d %d %d\n", i, j, k, l);
	double save_a = A(i,j,k,l);
	A(i,j,k,l) = change_value;

		for (int ii = min_i ; ii <= max_i ; ++ii)
		{
		for (int jj = min_j ; jj <= max_j ; ++jj)
		{
		for (int kk = min_k ; kk <= max_k ; ++kk)
		{
		for (int ll = min_l ; ll <= max_l ; ++ll)
		{
		if (print_flag)
		   then printf("      checking %d %d %d %d\n", ii, jj, kk, ll);

		double should_be = ((ii == i) && (jj == j)
				     && (kk == k) && (ll == l))
				   ? change_value : VALUE_4D(ii,jj,kk,ll);
		assert( fuzzy<double>::EQ(A(ii,jj,kk,ll), should_be) );

		double diff = A(ii,jj,kk,ll) - should_be;
		sumsq_of_diff += diff*diff;
		}
		}
		}
		}
	A(i,j,k,l) = save_a;
	}
	}
	}
	}
double RMS_diff = sqrt(sumsq_of_diff / A.N_array());

printf("==> everything looks ok (RMS_diff=%.3g)\n", RMS_diff);
}