#include "CrystalFp.h"
#include "CrystalFpAnalysis.h"
#include "ElementSymbols.h"
#include "ReadPoscar.h"
#include "simpleopt/SimpleOpt.h"
#include <stdio.h>
#ifndef WIN32
#include <unistd.h>
#endif
#include <vector>
#include <string>
#include <stdlib.h>

void ShowUsage(void)
{
	fprintf(stderr, "Usage: CrystalFp [options] poscar_file [energy_file]\n");

	fprintf(stderr, "-t At1,At2,At3...            Atom types\n");					
	fprintf(stderr, "--elements At1,At2,At3...    Atom types\n");		
	fprintf(stderr, "-r {energy_delta}            Energy threshold from the minimum energy value\n");
	fprintf(stderr, "--threshold-from-min {energy_delta} Energy threshold from the minimum energy value\n");
	fprintf(stderr, "-e {max_energy}              Energy threshold\n");			
	fprintf(stderr, "--energy-threshold {max_energy} Energy threshold\n");
	fprintf(stderr, "-b {bin_size}                Diffraction bin size (default: 0.05)\n");				
	fprintf(stderr, "--bin-size {bin_size}        Diffraction bin size (default: 0.05)\n");		
	fprintf(stderr, "-p {peak_size}               Peak smoothing size (default: 0.02)\n");					
	fprintf(stderr, "--peak-size {peak_size}      Peak smoothing size (default: 0.02)\n");		
	fprintf(stderr, "-g {threshold}               Distance threshold for grouping (default: 0.01)\n");
	fprintf(stderr, "--grouping-threshold {threshold} Distance threshold for grouping (default: 0.01)\n");
	fprintf(stderr, "-k {k_value}                 Number of shared nearest neighbors for grouping\n");
	fprintf(stderr, "--k {k_value}                Number of shared nearest neighbors for grouping\n");
	fprintf(stderr, "-s {max_step}                Max number of steps to read\n");
	fprintf(stderr, "--max-step {max_step}        Max number of steps to read\n");
	fprintf(stderr, "--end-step {last_step}       Last step to read\n");
	fprintf(stderr, "--start-step {first_step}    First step to read\n");
	fprintf(stderr, "-c {cutoff_distance}         Force cutoff distance\n");
	fprintf(stderr, "--cutoff-distance {cutoff_distance} Force cutoff distance\n");
	fprintf(stderr, "--triangle-fix               Apply the triangular inequality fix\n"};		
	fprintf(stderr, "--compute-q4                 Compute the cristallinity factor Q4\n");			
	fprintf(stderr, "--compute-q6                 Compute the cristallinity factor Q6\n");	
	fprintf(stderr, "--compute-w4                 Compute the cristallinity factor W4\n");	
	fprintf(stderr, "--compute-w6                 Compute the cristallinity factor W6\n");	
	fprintf(stderr, "-?                           This help\n");			
	fprintf(stderr, "--help                       This help\n");
	fprintf(stderr, "--debug-level {level}		  Set the debug level\n");
	fprintf(stderr, "--fingerprint-method {idx}   Set the fingerprint method (position in the enum)\n");
	fprintf(stderr, "--distance-method {idx}      Set the distance method (position in the enum)\n");
	fprintf(stderr, "--classification-method {idx} Set the distance method (position in the enum)\n");
}

// Decode SimpleOpt error
static const char *GetLastErrorText(CSimpleOpt& args)
{
    switch(args.LastError())
	{
    case SO_SUCCESS:            return "Success";
    case SO_OPT_INVALID:        return "Unrecognized option";
    case SO_OPT_MULTIPLE:       return "Option matched multiple strings";
    case SO_ARG_INVALID:        return "Option does not accept argument";
    case SO_ARG_INVALID_TYPE:   return "Invalid argument format";
    case SO_ARG_MISSING:        return "Required argument is missing";
    case SO_ARG_INVALID_DATA:   return "Invalid argument data";
    default:                    return "Unknown error";
    }
}

int main(int argc, char *argv[])
{
	// Setup the command line parsing
	enum {
		OPT_TYPES,
		OPT_REV_THRESH,
		OPT_THRESH,
		OPT_BIN_SIZE,
		OPT_PEAK_SIZE,
		OPT_GROUPING_THRESH,
		OPT_K_VALUE,
		OPT_MAX_STEPS,
		OPT_MIN_STEP,
		OPT_CUTOFF_DIST,
		OPT_TRI_FIX,
		OPT_HELP,
		OPT_DEBUG,
		OPT_FP_METHOD,
		OPT_DIST_METHOD,
		OPT_CLASS_METHOD,
		OPT_COMPUTE_Q4,
		OPT_COMPUTE_Q6,
		OPT_COMPUTE_W4,
		OPT_COMPUTE_W6
	};

	CSimpleOpt::SOption g_rgOptions[] = {
		{ OPT_TYPES,			"-t",					   SO_REQ_SEP },
		{ OPT_TYPES,			"--elements",			   SO_REQ_SEP },
		{ OPT_REV_THRESH,		"-r",					   SO_REQ_SEP },
		{ OPT_REV_THRESH,		"--threshold-from-min",	   SO_REQ_SEP },
		{ OPT_THRESH,			"-e",					   SO_REQ_SEP },
		{ OPT_THRESH,			"--energy-threshold",	   SO_REQ_SEP },
		{ OPT_BIN_SIZE,			"-b",					   SO_REQ_SEP },
		{ OPT_BIN_SIZE,			"--bin-size",			   SO_REQ_SEP },
		{ OPT_PEAK_SIZE,		"-p",					   SO_REQ_SEP },
		{ OPT_PEAK_SIZE,		"--peak-size",			   SO_REQ_SEP },
		{ OPT_GROUPING_THRESH,	"-g",					   SO_REQ_SEP },
		{ OPT_GROUPING_THRESH,	"--grouping-threshold",    SO_REQ_SEP },
		{ OPT_K_VALUE,			"-k",					   SO_REQ_SEP },
		{ OPT_K_VALUE,			"--k",					   SO_REQ_SEP },
		{ OPT_MAX_STEPS,		"-s",					   SO_REQ_SEP },
		{ OPT_MAX_STEPS,		"--max-step",			   SO_REQ_SEP },
		{ OPT_MAX_STEPS,		"--end-step",			   SO_REQ_SEP },
		{ OPT_MIN_STEP,			"--start-step",			   SO_REQ_SEP },
		{ OPT_CUTOFF_DIST,		"-c",					   SO_REQ_SEP },
		{ OPT_CUTOFF_DIST,		"--cutoff-distance",	   SO_REQ_SEP },
		{ OPT_TRI_FIX,			"--triangle-fix",		   SO_NONE    },
		{ OPT_COMPUTE_Q4,		"--compute-q4",			   SO_NONE    },
		{ OPT_COMPUTE_Q6,		"--compute-q6",			   SO_NONE    },
		{ OPT_COMPUTE_W4,		"--compute-w4",			   SO_NONE    },
		{ OPT_COMPUTE_W6,		"--compute-w6",			   SO_NONE    },
		{ OPT_HELP,				"-?",					   SO_NONE    },
		{ OPT_HELP,				"--help",				   SO_NONE    },
		{ OPT_DEBUG,			"--debug-level",		   SO_REQ_SEP },
		{ OPT_FP_METHOD,		"--fingerprint-method",    SO_REQ_SEP },
		{ OPT_DIST_METHOD,		"--distance-method",	   SO_REQ_SEP },
		{ OPT_CLASS_METHOD,		"--classification-method", SO_REQ_SEP },
		SO_END_OF_OPTIONS
	};

    // Declare our options parser, pass in the arguments from main as well as our array of valid options.
    CSimpleOpt args(argc, argv, g_rgOptions);

	// All the values gathered from the command line options
	std::string	atom_types;
	int			max_step = 0;
	int			min_step = 1;
	float		energy_threshold;
	bool		has_reverse_energy_threshold = false;
	bool		has_energy_threshold = false;
	float		cutoff_distance;
	bool		has_cutoff_distance = false;
	float		diffr_bin_size = 0.05F;
	float		diffr_peak_size = 0.02F;
	int			K = 0;
	float		max_distance_for_grouping = 0.01F;
	bool		do_triangle_fix = false;
	int			debug_level = 0;
	FingerprintType fp_method    = NORMALIZED_DIFFRACTION_HISTOGRAM;
	MeasureType 	dist_method  = COSINE_DISTANCE;
	ClassifyType 	class_method = PSEUDO_SNN_GROUPING;
	bool		compute_q4 = false;
	bool		compute_q6 = false;
	bool		compute_w4 = false;
	bool		compute_w6 = false;

    // While there are arguments left to process
    while(args.Next())
	{
        if(args.LastError() != SO_SUCCESS)
		{
            fprintf(stderr, "Error: %s for: %s\n", GetLastErrorText(args), args.OptionText());
            return 1;
        }

		if(debug_level == 1)
		{
			fprintf(stderr, "Option, ID: %3d, Text: '%s', Argument: '%s'\n",
					args.OptionId(), args.OptionText(),
					args.OptionArg() ? args.OptionArg() : "");
		}

		switch(args.OptionId())
		{
		case OPT_TYPES:
			atom_types = args.OptionArg();
			break;

		case OPT_REV_THRESH:
			energy_threshold = (float)atof(args.OptionArg());
			has_reverse_energy_threshold = true;
			break;

		case OPT_THRESH:
			energy_threshold = (float)atof(args.OptionArg());
			has_energy_threshold = true;
			break;

		case OPT_BIN_SIZE:
			diffr_bin_size = (float)atof(args.OptionArg());
			break;

		case OPT_PEAK_SIZE:
			diffr_peak_size = (float)atof(args.OptionArg());
			break;

		case OPT_GROUPING_THRESH:
			max_distance_for_grouping = (float)atof(args.OptionArg());
			break;

		case OPT_K_VALUE:
			K = (int)atol(args.OptionArg());
			break;

		case OPT_MAX_STEPS:
			max_step = atoi(args.OptionArg());
			break;

		case OPT_MIN_STEP:
			min_step = atoi(args.OptionArg());
			break;

		case OPT_CUTOFF_DIST:
			cutoff_distance = (float)atof(args.OptionArg());
			has_cutoff_distance = true;
			break;

		case OPT_TRI_FIX:
			do_triangle_fix = true;
			break;

		case OPT_HELP:
			ShowUsage();
			return 0;

		case OPT_DEBUG:
			debug_level = atoi(args.OptionArg());
			break;

		case OPT_FP_METHOD:
			fp_method = static_cast<FingerprintType>(atoi(args.OptionArg()));
			break;

		case OPT_DIST_METHOD:
			dist_method = static_cast<MeasureType>(atoi(args.OptionArg()));
			break;

		case OPT_CLASS_METHOD:
			class_method = static_cast<ClassifyType>(atoi(args.OptionArg()));
			break;

		case OPT_COMPUTE_Q4:
			compute_q4 = true;
			break;

		case OPT_COMPUTE_Q6:
			compute_q6 = true;
			break;

		case OPT_COMPUTE_W4:
			compute_w4 = true;
			break;

		case OPT_COMPUTE_W6:
			compute_w6 = true;
			break;

		case OPT_PEARSON:
			compute_pearson = true;

		case OPT_DIST_DISTRIB:
			output_distance_distribution = true;
			file_distance_distribution = args.OptionArg();
			break;
		}
	}

	// Set filenames (Poscar file and energy file)
	const char* fnp = 0;
	const char* fne = 0;

	switch(args.FileCount())
	{
	case 0:
		ShowUsage();
		return 0;

	default:
		fnp = args.File(0);
		fne = args.File(1);
		break;

	case 1:
		fnp = args.File(0);
		break;
	}

	// Prepare the atom name list
	std::vector<int> atom_z;
	if(atom_types.size() > 0)
	{
		const char* separator = " ,;\t\n";
		std::string::size_type start;
		std::string::size_type end = 0;

		std::vector<std::string> names;
		for(;;)
		{
			start = atom_types.find_first_not_of(separator, end);
			if(start == std::string::npos) break;
			end = atom_types.find_first_of(separator, start);

			names.push_back(std::string(atom_types, start, end - start));
		}

		for(std::vector<std::string>::const_iterator ia = names.begin(); ia != names.end(); ++ia)
		{
			atom_z.push_back(ElementSymbolToZ(ia->c_str()));
		}
	}

	// For testing output the values read
	fprintf(stderr, "\n\n");
	if(fnp) fprintf(stderr, "POSCAR file: %s\n", fnp);
 	if(fne) fprintf(stderr, "ENERGY file: %s\n", fne);
	fprintf(stderr, "\n");

	if(atom_types.size() > 0)
	{
		fprintf(stderr, "Atoms types:");
		for(std::vector<int>::const_iterator ia = atom_z.begin(); ia != atom_z.end(); ++ia)
		{
			fprintf(stderr, " %s", ElementZToSymbol(*ia));
		}
		fprintf(stderr, "\n");
	}

	if(max_step > 0) fprintf(stderr, "Max step %d\n", max_step);
	fprintf(stderr, "Diffr. bin size   %.4f\n", diffr_bin_size);
	fprintf(stderr, "Diffr. peak size  %.4f\n", diffr_peak_size);
	fprintf(stderr, "K                 %d\n", K);
	fprintf(stderr, "Grouping distance %.4f\n", max_distance_for_grouping);
	if(do_step_fusion_progressing) fprintf(stderr, "Fusion progressing step %.4f\n", step_fusion_progressing);
	fprintf(stderr, "\n");

	// Debug of the argument processing. End here
	if(debug_level == 1) return 0;

	// Start the CrystalFp library
	CrystalFp sl;

	// Read requested steps
	if(ReadPoscarAndEnergies(fnp, fne, min_step, max_step, atom_z, sl)) return 1;

	// Select  based on energies
	if(sl.HasEnergies())
	{
		fprintf(stderr, "Min energy %.6f\n", sl.GetMinEnergy());

		// Select  based on energies
		if(has_energy_threshold)
		{
			fprintf(stderr, "Threshold energy %.6f\n", energy_threshold);
			sl.EnergyThreshold(energy_threshold);
		}
		else if(has_reverse_energy_threshold)
		{
			energy_threshold += sl.GetMinEnergy();
			fprintf(stderr, "Threshold energy %.6f\n", energy_threshold);
			sl.EnergyThreshold(energy_threshold);
		}
	}

	// Show how many structures has been loaded and selected
	fprintf(stderr, "Loaded %d of %d\n", sl.GetNumActiveStructures(), sl.GetNumTotalStructures());

	// Compute the number of NN to use
	if(has_cutoff_distance)
	{
		fprintf(stderr, "Forced cutoff distance: %.4f\n", cutoff_distance);
	}
	else
	{
		cutoff_distance = sl.ComputeCutoffDistance();
		fprintf(stderr, "Computed cutoff distance: %.4f\n", cutoff_distance);
	}

	// Compute the fingerprint
	int safe_nn = sl.ComputeFingerprints(fp_method, cutoff_distance, diffr_bin_size, diffr_peak_size);
	fprintf(stderr, "Fingerprint method: %s\n", sl.GetFingerprintMethod());
	fprintf(stderr, "Dimension: %d\n", safe_nn);

	// Compute the distance matrix
	sl.ComputeDistanceMatrix(dist_method);
	fprintf(stderr, "Computed distances using method: %s\n", sl.GetDistanceMethod());

	// Test triangle inequality
	CrystalFpAnalysis cfa(&sl);
	float violations = cfa.TriangleInequalityViolations();
	fprintf(stderr, "Triangle inequality violated: %.2f%%\n", 100.0F*violations);

	// Fix triangle inequalities
	if(do_triangle_fix)
	{
		fprintf(stderr, "\nFixing triangle inequality\n");
		int	sts = sl.FixTriangleInequalityViolations();
		if(sts) fprintf(stderr, "Status: %d\n", sts);
	}

	// Classify phase
	sl.ClassifyResults(class_method, max_distance_for_grouping, K);
	fprintf(stderr, "Grouping using method: %s\n", sl.GetGroupingMethod());
	fprintf(stderr, "N. groups %d - n. single %d\n", sl.GetNgroups(), sl.GetNsingle());

	// Compute the crystallinity factors Q4, Q6, W4 and W6
	if(compute_q4 || compute_q6 || compute_w4 || compute_w6)
	{
		int subparam;
		if(compute_q6)      subparam = 1;
		else if(compute_w4) subparam = 2;
		else if(compute_w6) subparam = 3;
		else                subparam = 0;

		CrystalFpAnalysis cfa(&sl, DEGREE_OF_CRYSTALLINITY, subparam);
		cfa.SetupAnalysis();

		int nval = cfa.GetNumValues();
		float *x = new float[nval];
		float *y = new float[nval];
		cfa.GetChartValues(x, y);

		for(int i=0; i < sl.GetNumActiveStructures(); ++i)
		{
			printf("%5.0f %16.11f\n", x[i], y[i]);
		}

		delete [] x;
		delete [] y;
	}
}