// test_fd_grid.cc -- test driver for fd_grid and lower-level classes
// $Id$

//
// main
/// test_grid - test various grid:: member functions
///
/// f - test function
/// derivs - analytical derivatives of f (code generated by Maple)
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

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

#include "fp.hh"
#include "coords.hh"
#include "grid.hh"
#include "fd_grid.hh"

using jtutil::error_exit;

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

//
// This program is a test driver for the classes fd_grid::, grid_arrays::,
// and grid::.
//

const char *usage =
"\
Usage:\n\
   test_fd_grid  delta_drho  delta_dsigma  [ N_rhosigma ]\n\
                 [ fn | drho | dsigma | drhorho | drhosigma | dsigmasigma ] \n\
where\n\
   delta_drho  and  delta_dsigma  are the (rho,sigma) grid spacings\n\
   N_rhosigma , if specified, shrinks the max_{rho,sigma} grid bounds so\n\
      the grid has about that many points in each of the rho and sigma\n\
      directions\n\
\n\
This program tests the fd_grid:: finite differencing code by setting\n\
up a test function, finite differencing it, computing the analytical\n\
derivatives, and printing a table of the difference between the finite\n\
difference and analytical derivatives.  It also tests various grid::\n\
member functions by comparing their results against correct results.\n\
\n\
By default, the finite differencing test uses a combination of all\n\
the partial derivatives.  If the 4th argument is specified, then this\n\
restricts the test to only that derivative (this is useful for debugging\n\
problems with class fd_grid::).\n\
\n\
The program prints an ASCII data file on standard output, suitable for\n\
the gnuplot `splot' command, giving the error\n\
   finite differencing - analytical\n\
";

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

// bit-mask flags to encode which derivative to test
static const int flag_fn          = 0x01,
		 flag_drho        = 0x02,
		 flag_dsigma      = 0x04,
		 flag_drhorho     = 0x08,
		 flag_drhosigma   = 0x10,
		 flag_dsigmasigma = 0x20,
		 flag_all         = 0x3f;

// when testing multiple derivatives, we combine them together
// with these weights (chosen from digits of pi)
static const fp weight_fn          = 3.1,
		weight_drho        = 4.1,
		weight_dsigma      = 5.9,
		weight_drhorho     = 2.6,
		weight_drhosigma   = 5.3,
		weight_dsigmasigma = 5.8;

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

// fixed parameters for the grid

const fp min_drho   = - 70.0;
      fp max_drho   = +110.0;		// may be shrunk (via N_rhosigma)
const fp min_dsigma = -100.0;
      fp max_dsigma = +120.0;		// may be shruna (via N_rhosigma)

const int   min_rho_N_ghost_zones = 3;
const int   max_rho_N_ghost_zones = 4;
const int min_sigma_N_ghost_zones = 2;
const int max_sigma_N_ghost_zones = 5;

const int N_gridfns = 4;

// gfns
const int gfn_fn = 0;
const int gfn_deriv_true = 1;
const int gfn_deriv_fd = 2;
const int gfn_error = 3;

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

//
// function prototypes
//

namespace {
void test_grid(const grid& g,
	       double delta_drho, double delta_dsigma,
	       int min_irho, int max_irho, int min_isigma, int max_isigma);

void setup_fn(grid& g);
void setup_deriv_true(grid& g, int flags);
void setup_deriv_fd(grid& g, int flags);

fp fn(fp rho, fp sigma);
fp deriv(fp rho, fp sigma, int flags);
	  };

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

int main(int argc, const char *argv[])
{
//
// ***** parse the arguments
//

fp delta_drho, delta_dsigma;
if (! (    ((argc == 3) || (argc == 4) || (argc == 5))
	&& (sscanf(argv[1], FP_SCANF_FORMAT, &delta_drho) == 1)
	&& (sscanf(argv[2], FP_SCANF_FORMAT, &delta_dsigma) == 1)    ) )
   then error_exit(ERROR_EXIT, "%s", usage);			/*NOTREACHED*/

if (argc >= 4)
   then {
	int N_rhosigma;
	if (! (sscanf(argv[3], "%d", &N_rhosigma) == 1) )
	   then error_exit(ERROR_EXIT, "%s", usage);		/*NOTREACHED*/

	// size specified ==> shrink the grid to that size
	max_drho = min_drho + fp(N_rhosigma)*delta_drho;
	max_dsigma = min_dsigma + fp(N_rhosigma)*delta_dsigma;
	}

int flags = flag_all;
if (argc >= 5)
   then {
	if      (STRING_EQUAL(argv[4], "fn")) then flags = flag_fn;
	else if (STRING_EQUAL(argv[4], "drho")) then flags = flag_drho;
	else if (STRING_EQUAL(argv[4], "dsigma")) then flags = flag_dsigma;
	else if (STRING_EQUAL(argv[4], "drhorho")) then flags = flag_drhorho;
	else if (STRING_EQUAL(argv[4], "drhosigma"))
	   then flags = flag_drhosigma;
	else if (STRING_EQUAL(argv[4], "dsigmasigma"))
	   then flags = flag_dsigmasigma;
	else	error_exit(ERROR_EXIT, "%s", usage);		/*NOTREACHED*/
	}

if (! (    fuzzy<fp>::is_integer(min_drho  /delta_drho  )
	&& fuzzy<fp>::is_integer(max_drho  /delta_drho  )
	&& fuzzy<fp>::is_integer(min_dsigma/delta_dsigma)
	&& fuzzy<fp>::is_integer(max_dsigma/delta_dsigma)    ) )
   then error_exit(ERROR_EXIT,
		   "grid spacings must divide grid min/max bounds!\n");
								/*NOTREACHED*/
const int min_irho   = round<fp>::to_integer(min_drho  /delta_drho  );
const int max_irho   = round<fp>::to_integer(max_drho  /delta_drho  );
const int min_isigma = round<fp>::to_integer(min_dsigma/delta_dsigma);
const int max_isigma = round<fp>::to_integer(max_dsigma/delta_dsigma);

struct grid_arrays::array_pars grid_array_pars;
grid_array_pars.  min_irho =   min_irho;
grid_array_pars.  max_irho =   max_irho;
grid_array_pars.min_isigma = min_isigma;
grid_array_pars.max_isigma = max_isigma;
grid_array_pars.  min_rho_N_ghost_zones =   min_rho_N_ghost_zones;
grid_array_pars.  max_rho_N_ghost_zones =   max_rho_N_ghost_zones;
grid_array_pars.min_sigma_N_ghost_zones = min_sigma_N_ghost_zones;
grid_array_pars.max_sigma_N_ghost_zones = max_sigma_N_ghost_zones;

struct grid::grid_pars grid_pars;
grid_pars.    min_drho =     min_drho;
grid_pars.  delta_drho =   delta_drho;
grid_pars.    max_drho =     max_drho;
grid_pars.  min_dsigma =   min_dsigma;
grid_pars.delta_dsigma = delta_dsigma;
grid_pars.  max_dsigma =   max_dsigma;

printf("# ");
fd_grid g(grid_array_pars, grid_pars, N_gridfns, true);

test_grid(g,
	  delta_drho, delta_dsigma,
	  min_irho, max_irho, min_isigma, max_isigma);

setup_fn(g);
setup_deriv_true(g, flags);
setup_deriv_fd(g, flags);
}

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

//
// This function does some basic tests on grid:: and grid_arrays::
// member functions.
//
namespace {
void test_grid(const grid& g,
	       double delta_drho, double delta_dsigma,
	       int min_irho, int max_irho, int min_isigma, int max_isigma)
{
assert( g.N_gridfns() == N_gridfns );

assert( fuzzy<fp>::EQ(g.min_drho(), min_drho) );
assert( fuzzy<fp>::EQ(g.delta_drho(), delta_drho) );
assert( fuzzy<fp>::EQ(g.max_drho(), max_drho) );
assert( fuzzy<fp>::EQ(g.min_dsigma(), min_dsigma) );
assert( fuzzy<fp>::EQ(g.delta_dsigma(), delta_dsigma) );
assert( fuzzy<fp>::EQ(g.max_dsigma(), max_dsigma) );

assert( g.min_irho() == min_irho );
assert( g.max_irho() == max_irho );
assert( g.min_isigma() == min_isigma );
assert( g.max_isigma() == max_isigma );

assert( g.full_grid__min_irho() == min_irho - min_rho_N_ghost_zones );
assert( g.full_grid__max_irho() == max_irho + max_rho_N_ghost_zones );
assert( g.full_grid__min_isigma() == min_isigma - min_sigma_N_ghost_zones );
assert( g.full_grid__max_isigma() == max_isigma + max_sigma_N_ghost_zones );
}
	  }

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

//
// This function sets up the test function for the finite differencing
// tests.
//
namespace {
void setup_fn(grid& g)
{
printf("# setting up function\n");
	for (int irho = g.min_irho() ; irho <= g.max_irho() ; ++irho)
	{
	for (int isigma = g.min_isigma() ; isigma <= g.max_isigma() ; ++isigma)
	{
	const fp rho = g.rho_of_irho(irho);
	const fp sigma = g.sigma_of_isigma(isigma);
	g.gridfn(gfn_fn, irho,isigma) = fn(rho,sigma);
	}
	}
}
	  }

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

//
// This function sets up the analytic derivative expression for the
// finite differencing tests.
//
namespace {
void setup_deriv_true(grid& g, int flags)
{
printf("# setting up analytic derivatives\n");
	for (int irho = g.min_irho() ; irho <= g.max_irho() ; ++irho)
	{
	for (int isigma = g.min_isigma() ; isigma <= g.max_isigma() ; ++isigma)
	{
	const fp rho = g.rho_of_irho(irho);
	const fp sigma = g.sigma_of_isigma(isigma);
	g.gridfn(gfn_fn, irho,isigma) = deriv(rho,sigma, flags);
	}
	}
}
	  }

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

//
// This function sets up the analytic derivative expression for the
// finite differencing tests.
//
namespace {
void setup_deriv_fd(grid& g, int flags)
{
printf("# computing finite-difference derivatives\n");
	for (int irho = g.min_irho() ; irho <= g.max_irho() ; ++irho)
	{
	for (int isigma = g.min_isigma() ; isigma <= g.max_isigma() ; ++isigma)
	{
	const fp rho = g.rho_of_irho(irho);
	const fp sigma = g.sigma_of_isigma(isigma);
	//g.gridfn(gfn_fn, irho,isigma) = 
	}
	}
}
	  }

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

//
// this is the test function for our finite differencing tests
//
namespace {
fp fn(fp rho, fp sigma)
{
return exp(sin(rho) + tanh(sigma));
}
	  };

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

//
// This function computes the sum of our linear combination of
// various (analytical) derivatives of fn().  The derivatives were
// machine-generated via Maple's codegen[C]() function.
//
// "test_fd_grid.maple" is the Maple input
// "test_fd_grid.out" is the Maple input; code here is cut-n-pasted
//		       from the Maple codegen[C]() output there
//
// Arguments:
// flags = A bit-mask of flags specifying which terms to use in the
//	   linear combination.
//
namespace {
fp deriv(fp rho, fp sigma, int flags)
{
fp sum = 0.0;

// function itself
if (flags & flag_fn)
   then {
	fp temp = fn(rho,sigma);
	sum += weight_fn * temp;
	}

// 1st derivatives
if (flags & flag_drho)
   then {
	fp temp = cos(rho)*exp((sin(rho)*cosh(sigma)+sinh(sigma))
		  / cosh(sigma));
	sum += weight_drho * temp;
	}
if (flags & flag_dsigma)
   then {
	fp temp = exp((sin(rho)*cosh(sigma)+sinh(sigma))/cosh(sigma))
		  / pow(cosh(sigma),2.0);
	sum += weight_dsigma * temp;
	}

// 2nd derivatives
if (flags & flag_drhorho)
   then {
	fp temp = -exp((sin(rho)*cosh(sigma)+sinh(sigma))/cosh(sigma))
		   * (sin(rho)-pow(cos(rho),2.0));
	sum += weight_drhorho * temp;
	}
if (flags & flag_drhosigma)
   then {
	fp temp = cos(rho)*exp((sin(rho)*cosh(sigma)+sinh(sigma))/cosh(sigma))
		  / pow(cosh(sigma),2.0);
	sum += weight_drhosigma * temp;
	}
if (flags & flag_dsigmasigma)
   then {
	fp temp = -exp((sin(rho)*cosh(sigma)+sinh(sigma))/cosh(sigma))
		   * (2.0*sinh(sigma)*cosh(sigma)-1.0)
		   / pow(cosh(sigma),4.0);
	sum += weight_dsigmasigma * temp;
	}

return sum;
}
	  };