//
//  mapio.c
//  precisionControl
//
//  Created by Brian Schlining on 2007-08-31.
//  Copyright 2007 MBARI. All rights reserved.
//
// Documentation is formated for HeaderDoc or HeaderBrowser

#include "mapio.h"


//TODO this function fails to print which file or directory doesn't exist, making its error message near useless.
int check_error(int status, struct mapsrc* src) {
	const char* fname = "check_error";
	int err;
	
	if(status != NC_NOERR) {
		const char* status_tostring = nc_strerror(status);
		//char *status_tostring = nc_strerror(status);
		
		
		fprintf(stderr, "MAPIO::%s %s\n", fname, status_tostring);
		
		
		//free(status_tostring);
		// TODO how do we want to handle this error gracefully?
		src->status = src->status | MAPSRC_FILL_FAILURE;
		err = MAPIO_READERROR;
		//free(status_tostring);
	} else {
		err = MAPIO_OK;
	}
	return err;
}

void mapsrc_fill(const char* file, struct mapsrc* src) {

	const char* fname = "mapsrc_fill";
	int err, i;
	double range[2];
	double delta;
	
	// We don't refill existign strctures unless they've been free'd first
	if(src->x != NULL || src->y != NULL) {
		fprintf(stderr,
				"MAPIO::%s: An attempt was made to reallocate data to a mapsrc structure that already contains data. This may be a memory leak. You should call 'mapsrc_free' first",
				fname);
	}
	
	err = nc_open(file, NC_NOWRITE, &(src->ncid));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	/* get x variable */
	err = nc_inq_dimid(src->ncid, "x", &(src->xdimid));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	err = nc_inq_dimlen(src->ncid, src->xdimid, &(src->xdimlen));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	src->x = (double*) malloc(sizeof(double) * src->xdimlen);
	err = nc_inq_varid(src->ncid, "x", &(src->xid));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	err = nc_get_att_double(src->ncid, src->xid, "actual_range", range);
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	//fill in xvec based on range boundaries.
	delta = (range[1] - range[0]) / (src->xdimlen - 1);
	for(i = 0; i < int(src->xdimlen); i++) {
		src->x[i] = range[0] + delta * i;
	}
	
	//OLD WAY TO GET X-VAR: USES X VALUES FROM GRD FILE
	/*err = nc_get_var_double(src->ncid, src->xid, src->x);
	  if (MAPIO_OK != check_error(err, src)) { return; }
	*/
	
	/* get y variable */
	err = nc_inq_dimid(src->ncid, "y", &(src->ydimid));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	err = nc_inq_dimlen(src->ncid, src->ydimid, &(src->ydimlen));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	src->y = (double*) malloc(sizeof(double) * src->ydimlen);
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	err = nc_inq_varid(src->ncid, "y", &(src->yid));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	err = nc_get_att_double(src->ncid, src->yid, "actual_range", range);
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	//fill in xvec based on range boundaries.
	delta = (range[1] - range[0]) / (src->ydimlen - 1);
	for(i = 0; i < int(src->ydimlen); i++) {
		src->y[i] = range[0] + delta * i;
	}
	
	//OLD WAY TO GET Y-VAR: USES Y VALUES FROM GRD FILE
	/*err = nc_get_var_double(src->ncid, src->yid, src->y);
	  if (MAPIO_OK != check_error(err, src)) { return; }
	*/
	
	/* get z variable */
	err = nc_inq_varid(src->ncid, "z", &(src->zid));
	if(MAPIO_OK != check_error(err, src)) {
		return;
	}
	
	src->status = src->status | MAPSRC_IS_FILLED;
	
}

float mapsrc_find(struct mapsrc* src, double x, double y) {
	float* z;
	float z_out;
	int mapdata_code, XI = 1, YI = 0;
	
	size_t start[2];        // For NetCDF access -> {x0, y0}
	size_t count[2];        // For NetCDF access -> {xdimlen, ydimlen}
	
	// Make sure the data is within the map bounds
	struct mapbounds* bounds = mapbounds_init();
	mapbounds_fill1(src, bounds);
	mapdata_code = mapbounds_contains(bounds, x, y);
	free(bounds);
	
	if(mapdata_code != MAPBOUNDS_OUT_OF_BOUNDS) {
		// Get the nearest point in our map
		start[XI] = nearest(x, src->x, src->xdimlen);
		start[YI] = nearest(y, src->y, src->ydimlen);
		count[XI] = 1;
		count[YI] = 1;
		z = (float*) malloc(count[YI] * count[XI] * sizeof(float));
		check_error(nc_get_vara_float(src->ncid, src->zid, start, count, z), src);
		z_out = *z;
		free(z);
	} else {
		z_out = NAN;
	}
	
	return z_out;
}


struct mapsrc* mapsrc_init(void) {
	struct mapsrc* src = (struct mapsrc*) malloc(sizeof(struct mapsrc));
	src->ncid = 0;
	src->x = NULL;
	src->xid = 0;
	src->xdimid = 0;
	src->xdimlen = 0;
	src->y = NULL;
	src->yid = 0;
	src->ydimlen = 0;
	src->ydimid = 0;
	src->zid = 0;
	src->status = MAPSRC_IS_EMPTY;
	return src;
}


void mapsrc_free(struct mapsrc* src) {
	if(src->x != NULL) {
		free(src->x);
		src->x = NULL;
	}
	if(src->y != NULL) {
		free(src->y);
		src->y = NULL;
	}
	if(src != NULL) {
		free(src);
		src = NULL;
	}
}

char* mapsrc_tostring(struct mapsrc* src) {

	char* str = (char*) malloc(180 * sizeof(char));
	char buf[100];
	//char *buf = malloc(100 * sizeof(char));
	sprintf(str, "mapsrc {\n\tncid = %i\n", src->ncid);
	sprintf(buf, "\txid = %i\n", src->xid);
	strcat(str, buf);
	sprintf(buf, "\txdimid = %i\n", src->xdimid);
	strcat(str, buf);
	sprintf(buf, "\txdimlen = %d\n", int(src->xdimlen));
	strcat(str, buf);
	sprintf(buf, "\tyid = %i\n", src->yid);
	strcat(str, buf);
	sprintf(buf, "\tydimid = %i\n", src->ydimid);
	strcat(str, buf);
	sprintf(buf, "\tydimlen = %d\n", int(src->ydimlen));
	strcat(str, buf);
	sprintf(buf, "\tzid = %i\n", src->zid);
	strcat(str, buf);
	sprintf(buf, "\tstatus = %i\n", src->status);
	strcat(str, buf);
	strcat(str, "}");
	return str;
}

int mapdata_fill(struct mapsrc* src, struct mapdata* data, double x,
				 double y, double xwidth, double ywidth) {
				 
	const char* pname = "mapdata_fill";
	
	double xmin, xmax, ymin, ymax;
	// start = { x0, y0}, count = {xdimlen, ydimlen}
	size_t start[2];        // For NetCDF access -> {x0, y0}
	size_t count[2];        // For NetCDF access -> {xdimlen, ydimlen}
	//int i;
	int err, XI = 1, YI = 0;
	float* array;
	
	
	// Calculate the position of each corner of the submap
	xmin = x - xwidth / 2;
	xmax = x + xwidth / 2;
	ymin = y - ywidth / 2;
	ymax = y + ywidth / 2;
	//fprintf(stdout, "xmin: %f, xmax: %f, ymin: %f, ymax, %f", xmin, xmax, ymin, ymax);
	
	if(MAPIO_DEBUG) {
		fprintf(stdout, "MAPIO::%s: xmin: %f, xmax: %f, ymin: %f, ymax, %f", pname, xmin, xmax, ymin, ymax);
		fprintf(stdout, "MAPIO::%s: Looking for nearest values\n", pname);
	}
	// Get the indices into the map array for the corners of our submap
	start[XI] = nearest(xmin, src->x, src->xdimlen);
	start[YI] = nearest(ymin, src->y, src->ydimlen);
	count[XI] = nearest(xmax, src->x, src->xdimlen) - start[XI] + 1;
	count[YI] = nearest(ymax, src->y, src->ydimlen) - start[YI] + 1;
	
	
	// Assign values to 'data'
	if(MAPIO_DEBUG) {
		fprintf(stdout, "MAPIO::%s: Allocating x and y array memory\n", pname);
	}
	
	data->xdimlen = count[XI];
	data->ydimlen = count[YI];
	data->xpts = (double*) malloc(data->xdimlen * sizeof(double));
	memcpy(data->xpts, src->x + start[XI], count[XI] * sizeof(double));
	data->ypts = (double*) malloc(data->ydimlen * sizeof(double));
	memcpy(data->ypts, src->y + start[YI], count[YI] * sizeof(double));
	data->xcenter = (data->xpts[data->xdimlen - 1] + data->xpts[0]) / 2.0;
	data->ycenter = (data->ypts[data->ydimlen - 1] + data->ypts[0]) / 2.0;
	
	// Allocate memory for the arrays
	if(MAPIO_DEBUG) {
		fprintf(stdout, "MAPIO::%s: Allocating z array memory", pname);
	}
	array = (float*) malloc(count[YI] * count[XI] * sizeof(float));
	if(array == NULL) {
		fprintf(stderr, "MAPIO::%s: Out of memory. Failed to allocate memory for a mapdata structure\n", pname);
		data->status = data->status | MAPDATA_FILL_FAILURE;
	}
	
	//    if (MAPIO_DEBUG) {
	//        fprintf(stdout, "MAPIO::%s: Allocating z subarry memory", pname);
	//    }
	//    for (i = 0; i < count[YI]; i++) {
	//        *(array + i) = (float*) malloc(count[XI] * sizeof(float));
	//        if (*(array + i) = NULL) {
	//            fprintf(stderr, "MAPIO::%s: Out of memory. Failed to allocate memory for a mapdata structure\n", pname);
	//            // TODO OUT OF MEMORY!! handle error
	//            data->status = data->status | MAPDATA_FILL_FAILURE;
	//            return;
	//        }
	//    }
	
	
	// TODO Is this right...how the heck do I pass in an array of arrays to netcdf?
	data->z = array;
	
	// Extract the data from the netcdf source file.
	if(MAPIO_DEBUG) {
		fprintf(stdout, "MAPIO::%s: Reading z from netcdf", pname);
	}
	err = nc_get_vara_float(src->ncid, src->zid, start, count, (float*) data->z);
	check_error(err, src);
	data->status = MAPDATA_IS_FILLED;
	
	// Debug output used for compring results with matlab 'truth'
	if(MAPIO_DEBUG) {
		fprintf(stdout, "MAPIO::%s: ---- TEST REPORT ----\n", pname);
		
		fprintf(stdout, "Using data from %i\n", src->ncid);
		
		fprintf(stdout, "-- About X\n");
		fprintf(stdout, "X contains %i elements.\nMin = %f, Max = %f\n",
				int(src->xdimlen),
				src->x[0],
				src->x[int(src->xdimlen) - 1]);
		fprintf(stdout, "You specified: Center = %f, Width = %f\n", x, xwidth);
		fprintf(stdout, "Submap: x[%i] = %f to x[%i] = %f, %i elements\n",
				int(start[XI]),
				src->x[start[XI]],
				int(start[XI]) + int(data->xdimlen) - 1,
				src->x[start[XI] + int(data->xdimlen) - 1],
				int(data->xdimlen));
				
		fprintf(stdout, "-- About Y\n");
		fprintf(stdout, "Y contains %i elements.\nMin = %f, Max = %f\n",
				int(src->ydimlen),
				src->y[0],
				src->y[src->ydimlen - 1]);
		fprintf(stdout, "You specified: Center = %f, Width = %f\n", y, ywidth);
		fprintf(stdout, "Submap: y[%i] = %f to y[%i] = %f, %i elements\n",
				int(start[YI]),
				src->y[start[YI]],
				int(start[YI]) + int(data->ydimlen) - 1,
				src->y[start[YI] + int(data->ydimlen) - 1],
				int(data->ydimlen));
		fprintf(stdout, "\n");
		z_print((float*) data->z, data->ydimlen, data->xdimlen);
	}
	
	return mapdata_check(data, src, x, y, xwidth, ywidth);
}

void mapdata_free(struct mapdata* data, int free_all) {

	// Free the Z array
	//int i = 0;
	// TODO Confirm that I don' need to deallocate these memory blocks.
	//for (i = 0; i < data->xdimlen; i++) {
	//        printf("freeing z[%i] of %i\n", i, data->xdimlen);
	//        free(data->z[i]);
	//        data->z[i] = NULL;
	//    }
	if(data->z != NULL) {
		free(data->z);
		data->z = NULL;
	}
	
	if(data->xpts != NULL) {
		free(data->xpts);
		data->xpts = NULL;
	}
	
	if(data->ypts != NULL) {
		free(data->ypts);
		data->ypts = NULL;
	}
	
	if(free_all) {
		free(data);
	} else {
		data->xcenter = NAN;
		data->ycenter = NAN;
		data->xdimlen = 0;
		data->ydimlen = 0;
		data->status = MAPDATA_IS_EMPTY;
	}
	
}


char* mapdata_tostring(struct mapdata* data) {
	char* str = (char*) malloc(180 * sizeof(char));
	char buf[100];
	strcpy(str, "mapdata {\n");
	sprintf(buf, "\txcenter = %f\n", data->xcenter);
	strcat(str, buf);
	sprintf(buf, "\tycenter = %f\n", data->ycenter);
	strcat(str, buf);
	if(data->xpts == NULL) {
		sprintf(buf, "\tWARNING: x = NULL");
		strcat(str, buf);
	}
	if(data->ypts == NULL) {
		sprintf(buf, "\tWARNING: y = NULL");
		strcat(str, buf);
	}
	sprintf(buf, "\txdimlen = %u\n", int(data->xdimlen));
	strcat(str, buf);
	sprintf(buf, "\tydimlen = %u\n", int(data->ydimlen));
	strcat(str, buf);
	sprintf(buf, "\tstatus = %i\n", data->status);
	strcat(str, buf);
	strcat(str, "}");
	return str;
}

int nearest(double key, const double* base, size_t nmemb) {

	int i = 0, idx = 0, j;
	double a, dt, dt0;//, minval, maxval;
	//maxval = *(base + (nmemb - 1));
	//minval = *base;
	
	if(key > *(base + (nmemb - 1))) {
		idx = nmemb - 1;
		//fprintf(stderr, "MAPIO: Unable to find the nearest value to %f. It is more than the largest value, %f, in the array\n", key, maxval);
	} else if(key < *base) {
		idx = 0;
		//fprintf(stderr, "MAPIO: Unable to find the nearest value to %f. It is less than the smallest value, %f, in the array\n", key, minval);
	} else {
		dt0 = *(base + (nmemb - 1));
		for(j = 0; j < int(nmemb); j++) {
			a = *(base + j);
			dt = fabs(key - a);
			if(dt <= dt0) {
				dt0 = dt;
				idx = i;
			} else {
				break;
			}
			++i;
		}
	}
	return idx;
}

float getZ(const float* z, int row, int column, int columns) {
	return z[row * columns + column];
}

void z_print(const float* z, int rows, int columns) {
	float value;
	int i, j;
	
	for(i = 0; i < rows; i++) {
		for(j = 0; j < columns; j++) {
			//value = z[i * columns + j];
			value = getZ(z, i, j, columns);
			fprintf(stdout, "%f\t", value);
		}
		fprintf(stdout, "\n");
	}
	
}

int mapdata_check(struct mapdata* data, struct mapsrc* src, double xcenter, double ycenter, double xwidth, double ywidth) {

	int mapdata_code = MAPBOUNDS_OK;
	
	struct mapbounds* bounds = mapbounds_init();
	mapbounds_fill1(src, bounds);
	mapdata_code = mapbounds_contains(bounds, xcenter, ycenter);
	if(mapdata_code == MAPBOUNDS_OK) {
		mapdata_code = mapdata_checksize(bounds, data, xwidth, ywidth);
	}
	free(bounds);
	return mapdata_code;
}

struct mapbounds* mapbounds_init() {
	struct mapbounds* bounds = (struct mapbounds*) malloc(sizeof(struct mapbounds));
	bounds->dx = 0;
	bounds->dy = 0;
	bounds->ncid = 0;
	bounds->xmax = 0;
	bounds->xmin = 0;
	bounds->ymax = 0;
	bounds->ymin = 0;
	return bounds;
}


int mapbounds_fill1(struct mapsrc* src, struct mapbounds* bounds) {
	int return_code = MAPIO_OK;
	bounds->ncid = src->ncid;
	bounds->xmin = src->x[0];
	bounds->xmax = src->x[src->xdimlen - 1];
	bounds->dx = (bounds->xmax - bounds->xmin) / (float)(src->xdimlen - 1);
	bounds->ymin = src->y[0];
	bounds->ymax = src->y[src->ydimlen - 1];
	bounds->dy = (bounds->ymax - bounds->ymin) / (float)(src->ydimlen - 1);
	return return_code;
}

char* mapbounds_tostring(struct mapbounds* bounds) {
	char* str = (char*) malloc(256 * sizeof(char));
	char buf[100];
	strcpy(str, "mapbounds {\n");
	sprintf(buf, "\tncid = %d\n", bounds->ncid);
	strcat(str, buf);
	sprintf(buf, "\txmin = %f\n", bounds->xmin);
	strcat(str, buf);
	sprintf(buf, "\txmax = %f\n", bounds->xmax);
	strcat(str, buf);
	sprintf(buf, "\tdx = %f\n", bounds->dx);
	strcat(str, buf);
	sprintf(buf, "\tymin = %f\n", bounds->ymin);
	strcat(str, buf);
	sprintf(buf, "\tymax = %f\n", bounds->ymax);
	strcat(str, buf);
	sprintf(buf, "\tdy = %f\n", bounds->dy);
	strcat(str, buf);
	strcat(str, "}");
	return str;
}


int mapbounds_fill2(const char* file, struct mapbounds* bounds) {
	int return_code;
	struct mapsrc* src = mapsrc_init();
	mapsrc_fill(file, src);
	return_code = mapbounds_fill1(src, bounds);
	mapsrc_free(src);
	return return_code;
}

int mapbounds_contains(struct mapbounds* bounds, const double x, const double y) {
	int contains = MAPBOUNDS_OUT_OF_BOUNDS; // Default is False
	
	if(((double) bounds->xmax > x) &&
			((double) bounds->xmin < x) &&
			((double) bounds->ymax > y) &&
			((double) bounds->ymin < y)) {
		contains = MAPBOUNDS_OK;
	}
	
	return contains;
}

int mapdata_checksize(struct mapbounds* bounds, struct mapdata* data, const double xwidth, const double ywidth) {
	double x_size, y_size; // Number of expected pixels
	int size_is_ok = MAPBOUNDS_OK; // Default is True
	
	// Calculate expected size
	x_size = xwidth / bounds->dx;
	y_size = ywidth / bounds->dy;
	
	if(((double) data->xdimlen < (x_size - 1)) || ((double) data->ydimlen < (y_size - 1))) {
		size_is_ok = MAPBOUNDS_NEAR_EDGE;
	}
	return size_is_ok;
}