// io.cc -- I/O routines for this thorn
// $Header$
//
// input_gridfn - read an angular grid function from an input file
// output_gridfn - write an angular grid function to an output file
// output_Jacobians - write a Jacobian matrix or matrices to an output file
//
/// io_file_name - compute file name for angular-gridfn I/O file
//

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

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"

#include "../elliptic/Jacobian.hh"

#include "../gr/gfns.hh"
#include "../gr/gr.hh"

#include "horizon_sequence.hh"
#include "BH_diagnostics.hh"
#include "driver.hh"

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

//
// prototypes for functions local to this file
//

namespace {
const char* io_file_name(const struct IO_info& IO_info,
			 const char base_file_name[],
			 int hn, int AHF_iteration = 0);
	  }

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

//
// This function inputs a gridfn from a data file.
//
// We assume that this gridfn is ghosted, but the ghost zones are *not*
// present in the data file.
//
void input_gridfn(patch_system& ps, int unknown_gfn,
		  const struct IO_info& IO_info, const char base_file_name[],
		  int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
{
const char* file_name = io_file_name(IO_info, base_file_name,
				     hn, AHF_iteration);
if (print_msg_flag)
   then {
	if ( (unknown_gfn == gfns::gfn__h)
	     && (IO_info.time_iteration == 0) && (AHF_iteration == 0) )
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "   reading initial guess from \"%s\"", file_name);
	   else CCTK_VInfo(CCTK_THORNSTRING,
			   "   reading \"%s\"", file_name);
	}

switch	(IO_info.horizon_file_format)
	{
case horizon_file_format__ASCII_gnuplot:
	ps.read_ghosted_gridfn(unknown_gfn,
			       file_name,
			       false);		// no ghost zones in data file
	break;

case horizon_file_format__HDF5:
	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"input_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/

default:
	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   input_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n"
"                   (this should never happen!)"
		   ,
		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
	}
}

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

//
// This function outputs a gridfn from a data file.
//
// If the gridfn is h, then we also write out the xyz positions of the
// horizon surface points.
//
// FIXME: if the gridfn is not h, we assume that it's nominal-grid.
//
void output_gridfn(patch_system& ps, int unknown_gfn,
		   const struct IO_info& IO_info, const char base_file_name[],
		   int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
{
const char* file_name = io_file_name(IO_info, base_file_name,
				     hn, AHF_iteration);

switch	(IO_info.horizon_file_format)
	{
case horizon_file_format__ASCII_gnuplot:
	switch	(unknown_gfn)
		{
	case gfns::gfn__h:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "writing h to \"%s\"", file_name);
		ps.print_ghosted_gridfn_with_xyz
		      (unknown_gfn,
		       true, gfns::gfn__h,
		       file_name,
		       IO_info.output_ghost_zones_for_h); // should we include
							  // ghost zones?
		break;
	case gfns::gfn__Theta:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "writing Theta to \"%s\"", file_name);
		ps.print_gridfn(unknown_gfn, file_name);
		break;
	case gfns::gfn__Delta_h:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "writing Delta_h to \"%s\"", file_name);
		ps.print_gridfn(unknown_gfn, file_name);
		break;
	default:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "writing gfn=%d to \"%s\"",
				   unknown_gfn, file_name);
		ps.print_gridfn(unknown_gfn, file_name);
		break;
		}
	break;

case horizon_file_format__HDF5:
	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"output_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/

default:
	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   output_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n"
"                    (this should never happen!)"
		   ,
		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
	}
}

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

//
// This function prints one or two Jacobian matrices (and their difference
// in the latter case) to a named output file.
//
// Bugs:
// - The file format is hardwired to ASCII.
//
// Arguments:
// A null Jacobian pointer means to skip that Jacobian.
//
void output_Jacobians(const patch_system& ps,
                      const Jacobian* Jac_NP_ptr,
                      const Jacobian* Jac_SD_FDdr_ptr,
		      const struct IO_info& IO_info, const char base_file_name[],
		      int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
{
if (Jac_NP_ptr == NULL)
   then {
	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "output_Jacobians(): Jac_NP_ptr == NULL is not (yet) supported!");
	return;						// *** ERROR RETURN ***
	}

const char* file_name = io_file_name(IO_info, base_file_name,
				     hn, AHF_iteration);
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "   writing %s to \"%s\"",
		   ((Jac_SD_FDdr_ptr == NULL) ? "NP Jacobian"
					      : "NP and SD_FDdr Jacobians"),
		   file_name);

FILE *fileptr = fopen(file_name, "w");
if (fileptr == NULL)
   then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "output_Jacobians(): can't open output file \"%s\"!",
		   file_name);					/*NOTREACHED*/

fprintf(fileptr, "# column 1 = x patch name\n");
fprintf(fileptr, "# column 2 = x patch number\n");
fprintf(fileptr, "# column 3 = x irho\n");
fprintf(fileptr, "# column 4 = x isigma\n");
fprintf(fileptr, "# column 5 = x II\n");
fprintf(fileptr, "# column 6 = y patch name\n");
fprintf(fileptr, "# column 7 = y patch number\n");
fprintf(fileptr, "# column 8 = y irho\n");
fprintf(fileptr, "# column 9 = y isigma\n");
fprintf(fileptr, "# column 10 = y JJ\n");
fprintf(fileptr, "# column 11 = Jac_NP.element(II,JJ)\n");
if (Jac_SD_FDdr_ptr != NULL)
   then {
	fprintf(fileptr, "# column 12 = Jac_SD_FDdr.element(II,JJ)\n");
	fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n");
	fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n");
	}

    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
    {
    const patch& xp = ps.ith_patch(xpn);

    for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
    {
    for (int x_isigma = xp.min_isigma() ;
	 x_isigma <= xp.max_isigma() ;
	 ++x_isigma)
    {
    const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);

	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
	{
	const patch& yp = ps.ith_patch(ypn);

	for (int y_irho = yp.min_irho() ;
	     y_irho <= yp.max_irho() ;
	     ++y_irho)
	{
	for (int y_isigma = yp.min_isigma() ;
	     y_isigma <= yp.max_isigma() ;
	     ++y_isigma)
	{
	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);

	if (! Jac_NP_ptr->is_explicitly_stored(II,JJ))
	   then continue;			// skip sparse points

	const fp NP = Jac_NP_ptr->element(II,JJ);
	const fp SD_FDdr = (Jac_SD_FDdr_ptr == NULL)
			   ? 0.0
			   : Jac_SD_FDdr_ptr->element(II,JJ);

	if ((NP == 0.0) && (SD_FDdr == 0.0))
	   then continue;		// skip zero values (if == )

	fprintf(fileptr,
		"%s %d %d %d %d\t%s %d %d %d %d\t%#.10g",
		xp.name(), xpn, x_irho, x_isigma, II,
		yp.name(), ypn, y_irho, y_isigma, JJ,
		double(NP));

	if (Jac_SD_FDdr_ptr != NULL)
	   then {
		const fp abs_NP      = jtutil::abs(NP     );
		const fp abs_SD_FDdr = jtutil::abs(SD_FDdr);
		const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP);
		const fp SD_FDdr_abs_error = SD_FDdr - NP;
		const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;

		fprintf(fileptr,
			"\t%#.10g\t%e\t%e",
			double(SD_FDdr),
			double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
		}

	fprintf(fileptr, "\n");
	}
	}
	}
    }
    }
    }

fclose(fileptr);
}

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

//
// This function encapsulates our file-naming conventions for angular-gridfn
// output files (those used for h, H, and other angular grid functions).
//
// Arguments:
// base_file_name[] = from the parameter file
// hn = the horizon number
// AHF_iteration = the apparent horizon finder's internal iteration
//		   number (>= 1) if this is an intermediate iterate,
//		   or the default (0) if this is a final computed
//		   horizon position
//
// Results:
// This function returns (a pointer to) the file name.  The returned
// result points into a private static buffer; the usual caveats apply.
//
namespace {
const char* io_file_name(const struct IO_info& IO_info,
			 const char base_file_name[],
			 int hn, int AHF_iteration /* = 0 */)
{
const int N_file_name_buffer = 200;
static char file_name_buffer[N_file_name_buffer];

const char* file_name_extension;
switch	(IO_info.horizon_file_format)
	{
case horizon_file_format__ASCII_gnuplot:
	file_name_extension = IO_info.ASCII_gnuplot_file_name_extension;
	break;
case horizon_file_format__HDF5:
	file_name_extension = IO_info.HDF5_file_name_extension;
	break;
default:
	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   io_file_name(): unknown IO_info.horizon_file_format=(int)%d!\n"
"                   (this should never happen!)"
		   ,
		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
	}

if (AHF_iteration == 0)
   then snprintf(file_name_buffer, N_file_name_buffer,
		 "%s.t%d.ah%d.%s",
		 base_file_name,
		 IO_info.time_iteration, hn,
		 file_name_extension);
   else snprintf(file_name_buffer, N_file_name_buffer,
		 "%s.t%d.ah%d.it%d.%s",
		 base_file_name,
		 IO_info.time_iteration, hn, AHF_iteration,
		 file_name_extension);

return file_name_buffer;
}
	  }