#include "Utility.h"
#include "BaseGraphClass.h"
#include "GraphClass.h"
#include "NSPDK_FeatureGenerator.h"
#include "vectors.h"
#include "gzstream.h"
#include "OpenBabelConverter.h"

using namespace std;

const string PROG_CREDIT = "SVMSGDNSPDK Vers. 1.1 (19 Oct 2012)\nAuthor: Fabrizio Costa costa@informatik.uni-freiburg.de";
const string SEP = "----------------------------------------------------------------";
const string TAB = "    ";
const string CITATIONS = "For the graph kernels see Fabrizio Costa, Kurt De Grave, ''Fast Neighborhood Subgraph Pairwise Distance Kernel'', Proceedings of the 27th International Conference on Machine Learning (ICML-2010), Haifa, Israel, 2010..\nThe code for Stochastic Gradient Descent SVM is adapted from http://leon.bottou.org/projects/sgd. Léon Bottou and Yann LeCun, ''Large Scale Online Learning'', Advances in Neural Information Processing Systems 16, Edited by Sebastian Thrun, Lawrence Saul and Bernhard Schölkopf, MIT Press, Cambridge, MA, 2004.\nThe embedding method is adapted from L.Chen, A.Buja ''Local Multidimensional Scaling for Nonlinear Dimension Reduction, Graph Drawing, and Proximity Analysis'', Journal of the American Statistical Association, 2009.";

FlagsService& The_FlagsService = FlagsService::get_instance();

enum ActionType {
	TRAIN, TEST, PARAMETERS_OPTIMIZATION, CROSS_VALIDATION, CONFIDENCE, LEARNING_CURVE, TEST_PART, FEATURE, FEATURE_PART, MATRIX, EMBED
};

enum InputFileType {
	SPARSE_VECTOR, GRAPH, MOLECULAR_GRAPH, SEQUENCE, TREE
};

//------------------------------------------------------------------------------------------------------------------------
// Available losses
#define HINGELOSS 1
#define SMOOTHHINGELOSS 2
#define SQUAREDHINGELOSS 3
#define LOGLOSS 10
#define LOGLOSSMARGIN 11

// Select loss: NOTE: the selection done in the makefile
//#define LOSS LOGLOSS
//#define LOSS HINGELOSS
//#define LOSS SMOOTHHINGELOSS

// Zero when no bias
// One when bias term
#define BIAS 1

inline
double loss(double z) {
#if LOSS == LOGLOSS
	if (z > 18)
	return exp(-z);
	if (z < -18)
	return -z;
	return log(1+exp(-z));
#elif LOSS == LOGLOSSMARGIN
	if (z > 18)
	return exp(1-z);
	if (z < -18)
	return 1-z;
	return log(1+exp(1-z));
#elif LOSS == SMOOTHHINGELOSS
	if (z < 0)
	return 0.5 - z;
	if (z < 1)
	return 0.5 * (1-z) * (1-z);
	return 0;
#elif LOSS == SQUAREDHINGELOSS
	if (z < 1)
	return 0.5 * (1 - z) * (1 - z);
	return 0;
#elif LOSS == HINGELOSS
	if (z < 1) return 1 - z;
	return 0;
#else
	return 0;
#endif
}

inline
double dloss(double z) {
#if LOSS == LOGLOSS
	if (z > 18)
	return exp(-z);
	if (z < -18)
	return 1;
	return 1 / (exp(z) + 1);
#elif LOSS == LOGLOSSMARGIN
	if (z > 18)
	return exp(1-z);
	if (z < -18)
	return 1;
	return 1 / (exp(z-1) + 1);
#elif LOSS == SMOOTHHINGELOSS
	if (z < 0)
	return 1;
	if (z < 1)
	return 1-z;
	return 0;
#elif LOSS == SQUAREDHINGELOSS
	if (z < 1)
	return (1 - z);
	return 0;
#else
	if (z < 1) return 1;
	return 0;
#endif
}

//------------------------------------------------------------------------------------------------------------------------
enum OptionsType {
	FLAG, LIST, REAL, INTEGER, POSITIVE_INTEGER, STRING
};

class Parameters {
	class ParameterType {
	public:
		string mShort;
		string mLong;
		string mDescription;
		OptionsType mTypeCode;
		string mValue;
		bool mSet;
		vector<string> mValuesList;

	public:
		ParameterType() :
				mShort(""), mLong(""), mDescription(""), mTypeCode(STRING), mValue(""), mSet(false) {
		}

		void Parse(vector<string>& aParameterList) {
			string shortopt = "-" + mShort;
			string longopt = "--" + mLong;
			for (unsigned i = 0; i < aParameterList.size(); ++i) {
				if (aParameterList[i] == shortopt || aParameterList[i] == longopt) {
					mSet = true;
					switch (mTypeCode) {
					case FLAG:
						mValue = "1";
						break;
					case REAL: {
						if (i + 1 >= aParameterList.size()) throw range_error("ERROR: expected a value for option " + shortopt + "(" + longopt + ")");
						mValue = aParameterList[i + 1];
						double value = strtod(mValue.c_str(), NULL);
						if (!value) throw range_error("ERROR: Value for option " + shortopt + " (" + longopt + ") must be of type real");
					}
						break;
					case INTEGER: {
						if (i + 1 >= aParameterList.size()) throw range_error("ERROR: expected a value for option " + shortopt + "(" + longopt + ")");
						mValue = aParameterList[i + 1];
						int value = strtol(mValue.c_str(), NULL, 10);
						if (!value && mValue != "0") throw range_error("ERROR: Value for option " + shortopt + " (" + longopt + ") must be of type integer");
					}
						break;
					case POSITIVE_INTEGER: {
						if (i + 1 >= aParameterList.size()) throw range_error("ERROR: expected a value for option " + shortopt + "(" + longopt + ")");
						mValue = aParameterList[i + 1];
						unsigned value = strtoul(mValue.c_str(), NULL, 10);
						if (!value && mValue != "0") throw range_error("ERROR: Value for option " + shortopt + " (" + longopt + ") must be of type positive integer");
					}
						break;
					case STRING: {
						if (i + 1 >= aParameterList.size()) throw range_error("ERROR: expected a value for option " + shortopt + " (" + longopt + ")");
						mValue = aParameterList[i + 1];
					}
						break;
					case LIST: {
						if (i + 1 >= aParameterList.size()) throw range_error("ERROR: expected a value for option " + shortopt + " (" + longopt + ")");
						mValue = aParameterList[i + 1];
						bool test = false;
						for (unsigned i = 0; i < mValuesList.size(); ++i) {
							if (mValuesList[i] == mValue) {
								test = true;
								break;
							}
						}
						if (test == false) {
							string value_list = "";
							for (unsigned i = 0; i < mValuesList.size() - 1; ++i)
								value_list += mValuesList[i] + " | ";
							value_list += mValuesList.back();
							throw range_error("ERROR: Value for option " + shortopt + " (" + longopt + ") must be one of: " + value_list + " instead of: " + mValue);
						}
					}
						break;
					default:
						throw range_error("ERROR: Option " + shortopt + " (" + longopt + ") must be set");
					}
				}
			}
		}

		void OutputCompact(ostream& out) {
			if (mLong == "" && mShort == "") throw range_error("ERROR: Option does not have short nor long switch");
			if (mShort != "") out << "-" << mShort << " ";
			if (mLong != "") out << "--" << mLong << " ";
			if (mTypeCode == LIST) {
				string value_list = "";
				for (unsigned i = 0; i < mValuesList.size() - 1; ++i)
					value_list += mValuesList[i] + " | ";
				value_list += mValuesList.back();
				out << value_list << " ";
			}
			if (mTypeCode != FLAG) if (mValue != "") out << "[" << mValue << "] ";
			if (mDescription != "") out << TAB << mDescription << " ";
			out << endl;
		}

		void OutputExtended(ostream& out) {
			if (mLong == "" && mShort == "") throw range_error("ERROR: Option does not have short nor long switch");
			if (mShort != "") out << "-" << mShort << endl;
			if (mLong != "") out << "--" << mLong << endl;
			out << TAB << "Data type: ";
			switch (mTypeCode) {
			case FLAG:
				out << "flag" << endl;
				break;
			case REAL:
				out << "real number" << endl;
				break;
			case INTEGER:
				out << "integer number" << endl;
				break;
			case POSITIVE_INTEGER:
				out << "positive integer number" << endl;
				break;
			case STRING:
				out << "string" << endl;
				break;
			case LIST: {
				string value_list = "";
				for (unsigned i = 0; i < mValuesList.size() - 1; ++i)
					value_list += mValuesList[i] + " | ";
				value_list += mValuesList.back();
				out << "One of the following " << mValuesList.size() << " alternatives: " << value_list << endl;
			}
				break;
			default:
				throw range_error("ERROR4.1: Unrecognized parameter type code [" + mTypeCode);
			}
			if (mTypeCode != FLAG) {
				if (mValue != "") out << TAB << "Default: " << mValue << endl;
			}
			if (mDescription != "") out << TAB << "Description: " << mDescription << endl;
			out << endl;
		}
	}
	;

public:
	vector<ParameterType> mOptionList;
	string mAction;
	ActionType mActionCode;
	string mInputDataFileName;
	string mTargetFileName;
	string mModelFileName;
	string mFileType;
	string mOpenBabelFormat;
	InputFileType mFileTypeCode;
	bool mKernelNoNormalization;
	bool mMinKernel;
	unsigned mRadius;
	unsigned mDistance;
	unsigned mVertexDegreeThreshold;
	double mLambda;
	int mEpochs;
	unsigned mHashBitSize;
	unsigned mCrossValidationNumFolds;
	unsigned mLearningCurveNumPoints;
	unsigned mRandomSeed;
	string mKernelType;
	string mGraphType;
	unsigned mSemiSupervisedNumIterations;
	double mSemiSupervisedThreshold;
	bool mSemiSupervisedInduceOnlyPositive;
	bool mSemiSupervisedInduceOnlyNegative;
	string mSuffix;
	unsigned mSequenceDegree;
	unsigned mLMDSNumRandomRestarts;
	unsigned mLMDSNumIterations;
	unsigned mLMDSDimensionality;
	double mLMDSIterationEpsilon;
	unsigned mLMDSNeighborhoodSize;
	unsigned mLMDSNonNeighborhoodSize;
	unsigned mLMDSNeighborhoodSizeRange;
	double mLMDSTau;
	unsigned mLMDSTauExponentRange;
	bool mVerbose;
	bool mSequenceToken;
	bool mSequenceMultiLine;
	bool mSequencePairwiseInteraction;

public:
	Parameters() {
		SetupOptions();
	}

	void SetupOptions() {
		{
			ParameterType param;
			param.mShort = "h";
			param.mLong = "help";
			param.mDescription = "Prints compact help.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "H";
			param.mLong = "Help";
			param.mDescription = "Prints extended help.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "a";
			param.mLong = "action";
			param.mDescription = "";
			param.mTypeCode = LIST;
			param.mValue = "FEATURE";
			param.mValuesList.push_back("TRAIN");
			param.mValuesList.push_back("TEST");
			param.mValuesList.push_back("TEST_PART");
			param.mValuesList.push_back("CROSS_VALIDATION");
			param.mValuesList.push_back("CONFIDENCE");
			param.mValuesList.push_back("PARAMETERS_OPTIMIZATION");
			param.mValuesList.push_back("LEARNING_CURVE");
			param.mValuesList.push_back("FEATURE");
			param.mValuesList.push_back("FEATURE_PART");
			param.mValuesList.push_back("MATRIX");
			param.mValuesList.push_back("EMBED");
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "f";
			param.mLong = "file_type";
			param.mDescription = "";
			param.mTypeCode = LIST;
			param.mValue = "GRAPH";
			param.mValuesList.push_back("SPARSE_VECTOR");
			param.mValuesList.push_back("GRAPH");
			param.mValuesList.push_back("MOLECULAR_GRAPH");
			param.mValuesList.push_back("SEQUENCE");
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "k";
			param.mLong = "kernel_type";
			param.mDescription = "";
			param.mTypeCode = LIST;
			param.mValue = "NSPDK";
			param.mValuesList.push_back("NSPDK");
			param.mValuesList.push_back("WDK");
			param.mValuesList.push_back("PBK");
			param.mValuesList.push_back("USPK");
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "g";
			param.mLong = "graph_type";
			param.mDescription = "";
			param.mTypeCode = LIST;
			param.mValue = "UNDIRECTED";
			param.mValuesList.push_back("DIRECTED");
			param.mValuesList.push_back("UNDIRECTED");
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "no_normalization";
			param.mDescription = "Kernel parameter.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "min_kernel";
			param.mDescription = "Kernel parameter.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "o";
			param.mLong = "open_babel_file_format";
			param.mDescription = "MOLECULAR_GRAPH parameter.";
			param.mTypeCode = STRING;
			param.mValue = "sdf";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "i";
			param.mLong = "input_data_file_name";
			param.mDescription = "";
			param.mTypeCode = STRING;
			param.mValue = "";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "t";
			param.mLong = "target_file_name";
			param.mDescription = "";
			param.mTypeCode = STRING;
			param.mValue = "";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "m";
			param.mLong = "model_file_name";
			param.mDescription = "";
			param.mTypeCode = STRING;
			param.mValue = "model";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "s";
			param.mLong = "suffix";
			param.mDescription = "Suffix string for all output files.";
			param.mTypeCode = STRING;
			param.mValue = "";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "r";
			param.mLong = "radius";
			param.mDescription = "Kernel parameter.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "2";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "d";
			param.mLong = "distance";
			param.mDescription = "Kernel parameter.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "5";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "v";
			param.mLong = "vertex_degree_threshold";
			param.mDescription = "Kernel parameter.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "7";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "M";
			param.mLong = "sequence_degree";
			param.mDescription = "SEQUENCE data type.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "1";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "b";
			param.mLong = "hash_bit_size";
			param.mDescription = "Kernel parameter.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "14";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "l";
			param.mLong = "lambda";
			param.mDescription = "Stochastic gradient descend algorithm.";
			param.mTypeCode = REAL;
			param.mValue = "1e-6";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "e";
			param.mLong = "epochs";
			param.mDescription = "Stochastic gradient descend algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "5";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "c";
			param.mLong = "num_cross_validation_folds";
			param.mDescription = "";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "10";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "p";
			param.mLong = "num_learning_curve_evaluation_points";
			param.mDescription = "";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "10";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "R";
			param.mLong = "random_seed";
			param.mDescription = "";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "1";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "S";
			param.mLong = "num_iterations";
			param.mDescription = "In semi-supervised setting.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "3";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "T";
			param.mLong = "threshold";
			param.mDescription = "In semi-supervised setting. Only the top and low quantile will be used as positives and negative instances. A threshold of 1 means that all unsupervised instaces are used in the next phase.";
			param.mTypeCode = REAL;
			param.mValue = "1";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "only_positive";
			param.mDescription = "In semi-supervised setting. Induce only positive class instances.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "only_negative";
			param.mDescription = "In semi-supervised setting. Induce only negative class instances.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "N";
			param.mLong = "num_of_random_restarts";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "1";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "I";
			param.mLong = "num_of_iterations";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "5000";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "D";
			param.mLong = "dimensionality";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "2";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "E";
			param.mLong = "epsilon";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = REAL;
			param.mValue = "0.01";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "K";
			param.mLong = "neighborhood_size";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "10";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "n";
			param.mLong = "non_neighborhood_size";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "100";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "G";
			param.mLong = "neighborhood_size_range";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "U";
			param.mLong = "tau";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = REAL;
			param.mValue = "0.0005";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "u";
			param.mLong = "tau_exponent_range";
			param.mDescription = "In low dimensionality embedding algorithm.";
			param.mTypeCode = POSITIVE_INTEGER;
			param.mValue = "1";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "V";
			param.mLong = "verbose";
			param.mDescription = "Outputs the graphs and the feature encodings.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "sequence_token";
			param.mDescription = "Labels are strings separated by spaces.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "sequence_multi_line";
			param.mDescription = "The annotation is encoded on subsequent rows.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}
		{
			ParameterType param;
			param.mShort = "";
			param.mLong = "sequence_pairwise_interaction";
			param.mDescription = "Abstraction vertices relating all pairs of disjointed sequences are inserted.";
			param.mTypeCode = FLAG;
			param.mValue = "0";
			mOptionList.push_back(param);
		}

	}

	void Usage(string aCommandName, string aCompactOrExtended) {
		cerr << "PROGRAM: " << aCommandName << endl;
		cerr << "OPTIONS:" << endl;
		for (unsigned i = 0; i < mOptionList.size(); ++i)
			if (aCompactOrExtended == "EXTENDED") mOptionList[i].OutputExtended(cerr);
			else mOptionList[i].OutputCompact(cerr);
		cout << SEP << endl << CITATIONS << endl << SEP << endl;
		exit(0);
	}

	void Init(int argc, const char** argv) {
		if (argc == 1) {
			cerr << "Use -h for compact help and -H for extended help." << endl;
			exit(1);
		}

		//convert argc in an option string vector
		vector<string> options;
		for (int i = 1; i < argc; i++)
			options.push_back(argv[i]);

		//parse the option string vector
		for (unsigned i = 0; i < mOptionList.size(); ++i)
			mOptionList[i].Parse(options);

		//set the boolean parmeters to a default value of false
		mKernelNoNormalization = false;
		mMinKernel = false;
		mSemiSupervisedInduceOnlyPositive = false;
		mSemiSupervisedInduceOnlyNegative = false;
		mVerbose = false;
		mSequenceToken = false;
		mSequenceMultiLine = false;
		mSequencePairwiseInteraction = false;

		//set the data members of Parameters according to user choice
		for (unsigned i = 0; i < mOptionList.size(); ++i) {
			if (mOptionList[i].mSet) {
				if (mOptionList[i].mShort == "h") {
					Usage(argv[0], "COMPACT");
					exit(1);
				}
				if (mOptionList[i].mShort == "H") {
					Usage(argv[0], "EXTENDED");
					exit(1);
				}
				if (mOptionList[i].mShort == "i") mInputDataFileName = mOptionList[i].mValue;
				if (mOptionList[i].mShort == "t") mTargetFileName = mOptionList[i].mValue;
				if (mOptionList[i].mLong == "no_normalization") mKernelNoNormalization = true;
				if (mOptionList[i].mLong == "min_kernel") mMinKernel = true;
				if (mOptionList[i].mLong == "only_positive") mSemiSupervisedInduceOnlyPositive = true;
				if (mOptionList[i].mLong == "only_negative") mSemiSupervisedInduceOnlyNegative = true;
				if (mOptionList[i].mShort == "V") mVerbose = true;
				if (mOptionList[i].mLong == "sequence_token") mSequenceToken = true;
				if (mOptionList[i].mLong == "sequence_multi_line") mSequenceMultiLine = true;
				if (mOptionList[i].mLong == "sequence_pairwise_interaction") mSequencePairwiseInteraction = true;
			}
			if (mOptionList[i].mShort == "a") mAction = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "m") mModelFileName = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "f") mFileType = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "k") mKernelType = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "g") mGraphType = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "o") mOpenBabelFormat = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "s") mSuffix = mOptionList[i].mValue;
			if (mOptionList[i].mShort == "r") mRadius = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "d") mDistance = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "v") mVertexDegreeThreshold = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "M") mSequenceDegree = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "b") mHashBitSize = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "l") mLambda = stream_cast<double>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "e") mEpochs = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "c") mCrossValidationNumFolds = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "p") mLearningCurveNumPoints = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "R") mRandomSeed = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "S") mSemiSupervisedNumIterations = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "T") mSemiSupervisedThreshold = stream_cast<double>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "N") mLMDSNumRandomRestarts = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "I") mLMDSNumIterations = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "D") mLMDSDimensionality = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "E") mLMDSIterationEpsilon = stream_cast<double>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "K") mLMDSNeighborhoodSize = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "n") mLMDSNonNeighborhoodSize = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "G") mLMDSNeighborhoodSizeRange = stream_cast<unsigned>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "U") mLMDSTau = stream_cast<double>(mOptionList[i].mValue);
			if (mOptionList[i].mShort == "u") mLMDSTauExponentRange = stream_cast<unsigned>(mOptionList[i].mValue);
		}

		//check that set parameters are compatible
		if (mInputDataFileName == "") throw range_error("ERROR2.4: -i <input data file name> is missing.");

		if ((mAction == "TRAIN" || mAction == "CROSS_VALIDATION" || mAction == "LEARNING_CURVE") && mTargetFileName == "") throw range_error("ERROR2.5: -t <target file name> is missing.");

		//convert action string to action code
		if (mAction == "TRAIN") mActionCode = TRAIN;
		else if (mAction == "TEST") mActionCode = TEST;
		else if (mAction == "CROSS_VALIDATION") mActionCode = CROSS_VALIDATION;
		else if (mAction == "CONFIDENCE") mActionCode = CONFIDENCE;
		else if (mAction == "PARAMETERS_OPTIMIZATION") mActionCode = PARAMETERS_OPTIMIZATION;
		else if (mAction == "LEARNING_CURVE") mActionCode = LEARNING_CURVE;
		else if (mAction == "TEST_PART") mActionCode = TEST_PART;
		else if (mAction == "FEATURE") mActionCode = FEATURE;
		else if (mAction == "FEATURE_PART") mActionCode = FEATURE_PART;
		else if (mAction == "MATRIX") mActionCode = MATRIX;
		else if (mAction == "EMBED") mActionCode = EMBED;
		else throw range_error("ERROR2.46: Unrecognized action: <" + mAction + ">");

		//convert file type string to file type code
		if (mFileType == "GRAPH") mFileTypeCode = GRAPH;
		else if (mFileType == "SPARSE_VECTOR") mFileTypeCode = SPARSE_VECTOR;
		else if (mFileType == "MOLECULAR_GRAPH") mFileTypeCode = MOLECULAR_GRAPH;
		else if (mFileType == "SEQUENCE") mFileTypeCode = SEQUENCE;
		else if (mFileType == "TREE") mFileTypeCode = TREE;
		else throw range_error("ERROR2.46: Unrecognized file type: <" + mFileType + ">");
	}
};

//------------------------------------------------------------------------------------------------------------------------
class Kernel {
public:
	NSPDK_FeatureGenerator* mpFeatureGenerator;
	Parameters* mpParameters;
public:
	void Init(Parameters* apParameters) {
		mpParameters = apParameters;
		if (mpParameters->mKernelType == "NSPDK") mpFeatureGenerator = new ANSPDK_FeatureGenerator("anspdk");
		else if (mpParameters->mKernelType == "PBK") mpFeatureGenerator = new PBK_FeatureGenerator("pbk");
		else if (mpParameters->mKernelType == "WDK") mpFeatureGenerator = new WDK_FeatureGenerator("wdk");
		else if (mpParameters->mKernelType == "USPK") mpFeatureGenerator = new USPK_FeatureGenerator("uspk");
		else throw range_error("ERROR2.1: Unknown kernel type: " + mpParameters->mKernelType);

		ParametersSetup();
		cout << SEP << endl << "Kernel parameters" << endl << SEP << endl;
		mpFeatureGenerator->OutputParameters(cout);
		cout << SEP << endl;
	}

	void ParametersSetup() {
#ifdef DEBUGON
		mpFeatureGenerator->set_flag("verbosity", stream_cast<string>(1));
#endif
		if (mpParameters->mVerbose) mpFeatureGenerator->set_flag("verbosity", stream_cast<string>(1));
		if (mpParameters->mMinKernel) mpFeatureGenerator->set_flag("min_kernel", "true");
		if (mpParameters->mKernelNoNormalization) mpFeatureGenerator->set_flag("normalization", "false");
		mpFeatureGenerator->set_flag("radius", stream_cast<string>(mpParameters->mRadius));
		mpFeatureGenerator->set_flag("distance", stream_cast<string>(mpParameters->mDistance));
		mpFeatureGenerator->set_flag("hash_bit_size", stream_cast<string>(mpParameters->mHashBitSize));
		unsigned bitmask = (2 << mpParameters->mHashBitSize) - 1;
		mpFeatureGenerator->set_flag("hash_bit_mask", stream_cast<string>(bitmask));

		//if type of kernel PBK then also perform the following initializations
		if (mpParameters->mKernelType == "PBK") {
			mpFeatureGenerator->set_flag("lower_vertex_degree_threshold", stream_cast<string>(mpParameters->mVertexDegreeThreshold));
			mpFeatureGenerator->set_flag("vertex_degree_threshold", stream_cast<string>(mpParameters->mVertexDegreeThreshold));

			if (mpParameters->mMinKernel) dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("min_kernel", "true");
			if (mpParameters->mKernelNoNormalization) dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("normalization", "false");
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("radius", stream_cast<string>(mpParameters->mRadius));
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("distance", stream_cast<string>(mpParameters->mDistance));
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("hash_bit_size", stream_cast<string>(mpParameters->mHashBitSize));
			bitmask = (2 << mpParameters->mHashBitSize) - 1;
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("hash_bit_mask", stream_cast<string>(bitmask));
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mWDK.set_flag("lower_vertex_degree_threshold", stream_cast<string>(mpParameters->mVertexDegreeThreshold));

			if (mpParameters->mMinKernel) dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("min_kernel", "true");
			if (mpParameters->mKernelNoNormalization) dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("normalization", "false");
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("radius", stream_cast<string>(mpParameters->mRadius));
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("distance", stream_cast<string>(mpParameters->mDistance));
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("hash_bit_size", stream_cast<string>(mpParameters->mHashBitSize));
			bitmask = (2 << mpParameters->mHashBitSize) - 1;
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("hash_bit_mask", stream_cast<string>(bitmask));
			dynamic_cast<PBK_FeatureGenerator*>(mpFeatureGenerator)->mANSPDK.set_flag("vertex_degree_threshold", stream_cast<string>(mpParameters->mVertexDegreeThreshold));
		}
	}

	void GenerateFeatureVector(GraphClass& aG, SVector& oX) {
		mpFeatureGenerator->generate_feature_vector(aG, oX);
	}

	void GenerateVertexFeatureVector(GraphClass& aG, vector<SVector>& oXList) {
		mpFeatureGenerator->generate_vertex_feature_vector(aG, oXList);
	}

	double ComputeKernel(GraphClass& aG, GraphClass& aM) {
		SVector xg;
		SVector xm;
		GenerateFeatureVector(aG, xg);
		GenerateFeatureVector(aM, xm);
		return ComputeKernel(xg, xm);
	}

	double ComputeKernel(SVector& aX, SVector& aZ) {
		return dot(aX, aZ);
	}
};

//------------------------------------------------------------------------------------------------------------------------
class Data {
public:
	Parameters* mpParameters;
	Kernel mKernel;
	vector<double> mTargetList;
	vector<SVector> mVectorList;
public:
	Data() {
	}

	void Init(Parameters* apParameters) {
		mpParameters = apParameters;
		mKernel.Init(mpParameters);
	}

	double ComputeKernel(unsigned i, unsigned j) {
		if (i > Size() || j > Size()) throw range_error("ERROR3.1: Kernel computation for instances out of range. Available range: 0-" + stream_cast<string>(Size()) + " but asked for kernel between instances with index " + stream_cast<string>(i) + "," + stream_cast<string>(j));
		return mKernel.ComputeKernel(mVectorList[i], mVectorList[j]);
	}

	void LoadTarget() {
		mTargetList.clear();
		cout << "Reading target file: " << mpParameters->mTargetFileName << " ..";
		ifstream fin;
		fin.open(mpParameters->mTargetFileName.c_str());
		if (!fin) throw range_error("ERROR2.23: Cannot open file:" + mpParameters->mTargetFileName);
		while (!fin.eof()) {
			string line;
			getline(fin, line);
			stringstream ss;
			ss << line << endl;
			while (!ss.eof()) {
				double target(0);
				ss >> target;
				if (ss.good()) {
					assert(target==1 || target==-1 || target==0);
					mTargetList.push_back(target);
				}
			}
		}
		fin.close();
		cout << ".. read: " << mTargetList.size() << " targets." << endl;
	}

	void LoadGspanList(vector<GraphClass>& aGraphList) {
		mVectorList.clear();
		for (unsigned i = 0; i < aGraphList.size(); ++i) {
			SVector x;
			mKernel.GenerateFeatureVector(aGraphList[i], x);
			mVectorList.push_back(x);
		}
	}

	void LoadData() {
		mVectorList.clear();
		mKernel.ParametersSetup();
		igzstream fin;
		fin.open(mpParameters->mInputDataFileName.c_str());
		if (!fin) throw range_error("ERROR2.11: Cannot open file: " + mpParameters->mInputDataFileName);
		ProgressBar pb;
		cout << "Processing file: " << mpParameters->mInputDataFileName << endl;
		while (!fin.eof()) {
			switch (mpParameters->mFileTypeCode) {
			case GRAPH:
			case MOLECULAR_GRAPH:
			case SEQUENCE: {
				GraphClass g;
				SetGraphFromFile(fin, g);
				if (!g.IsEmpty()) {
					SVector x;
					mKernel.GenerateFeatureVector(g, x);
					mVectorList.push_back(x);
					pb.Count();
				}
			}
				break;
			case SPARSE_VECTOR: {
				SVector x;
				SetVectorFromSparseVectorAsciiFile(fin, x);
				if (x.size() > 0) {
					mVectorList.push_back(x);
					pb.Count();
				}
			}
				break;
			default:
				throw range_error("ERROR2.45: file type not recognized: " + mpParameters->mFileType);
			}
		}
	}

	void SetGraphFromFile(istream& in, GraphClass& oG) {
		switch (mpParameters->mFileTypeCode) {
		case GRAPH:
			SetGraphFromGraphGspanFile(in, oG);
			break;
		case MOLECULAR_GRAPH:
			SetGraphFromGraphOpenBabelFile(in, oG, mpParameters->mOpenBabelFormat);
			break;
		case SEQUENCE: {
			if (mpParameters->mSequenceToken && mpParameters->mSequenceMultiLine) SetGraphFromSequenceMultiLineTokenFile(in, oG);
			if (mpParameters->mSequenceToken && !mpParameters->mSequenceMultiLine) SetGraphFromSequenceTokenFile(in, oG);
			if (!mpParameters->mSequenceToken && mpParameters->mSequenceMultiLine) SetGraphFromSequenceMultiLineFile(in, oG);
			if (!mpParameters->mSequenceToken && !mpParameters->mSequenceMultiLine) SetGraphFromSequenceFile(in, oG);
			break;
		}
		default:
			throw range_error("ERROR2.45: file type not recognized: " + mpParameters->mFileType);
		}
	}

	void SetGraphFromGraphOpenBabelFile(istream& in, GraphClass& oG, string aFormat) {
		if (in.good() && !in.eof()) {
			OpenBabelConverter molecule_converter;
			molecule_converter.ConvertOpenBabelFormatToGraph(&in, oG, aFormat);
		}
	}

	void SetVectorFromSparseVectorAsciiFile(istream& in, SVector& aX) {
		string line;
		getline(in, line);
		if (line == "") return;
		stringstream ss;
		ss << line << endl;
		while (!ss.eof() && ss.good()) {
			string key_value;
			ss >> key_value;
			size_t limit = key_value.find_first_of(":", 0);
			if (limit != string::npos) { //if the delimiter ':' is found then proceed
				string key = key_value.substr(0, limit);
				string value = key_value.substr(limit + 1, key_value.size());
				unsigned key_int = stream_cast<unsigned>(key);
				double val_real = stream_cast<double>(value);
				aX.set(key_int, val_real);
			}
		}
	}

	void SetGraphFromSequenceFile(istream& in, GraphClass& oG) {
		vector<vector<unsigned> > vertex_component_list;
		vector<unsigned> vertex_component;
		bool graph_disconnect = true;
		string line;
		getline(in, line);
		if (line == "") return;
		unsigned vertex_counter = 0;
		for (unsigned position_counter = 0; position_counter < line.length(); position_counter++) {
			char char_label = line.at(position_counter);
			string label = stream_cast<string>(char_label);
			if (label == "|") {
				graph_disconnect = true;
				vertex_component_list.push_back(vertex_component);
				vertex_component.clear();
			} else {
				vertex_component.push_back(vertex_counter);
				AddVertexAndEdgesForSequence(oG, label, vertex_counter, graph_disconnect);
				vertex_counter++;
				graph_disconnect = false;
			}
		}
		if (mpParameters->mGraphType == "DIRECTED") AddReverseGraph(oG);
		vertex_component_list.push_back(vertex_component);
		if (vertex_component_list.size() > 1) AddAbstractConnections(oG, vertex_component_list);
	}

	void SetGraphFromSequenceTokenFile(istream& in, GraphClass& oG) {
		vector<vector<unsigned> > vertex_component_list;
		vector<unsigned> vertex_component;
		bool graph_disconnect = true;
		string line;
		getline(in, line);
		if (line == "") return;
		stringstream ss;
		ss << line << endl;
		unsigned vertex_counter = 0;
		while (!ss.eof() && ss.good()) {
			//add vertex
			string label;
			ss >> label;
			if (label != "") {
				if (label == "|") {
					graph_disconnect = true;
					vertex_component_list.push_back(vertex_component);
					vertex_component.clear();
				} else {
					vertex_component.push_back(vertex_counter);
					AddVertexAndEdgesForSequence(oG, label, vertex_counter, graph_disconnect);
					vertex_counter++;
					graph_disconnect = false;
				}
			}
		}
		if (mpParameters->mGraphType == "DIRECTED") AddReverseGraph(oG);
		vertex_component_list.push_back(vertex_component);
		if (vertex_component_list.size() > 1) AddAbstractConnections(oG, vertex_component_list);
	}

	void SetGraphFromSequenceMultiLineFile(istream& in, GraphClass& oG) {
		vector<vector<unsigned> > vertex_component_list;
		vector<unsigned> vertex_component;
		bool graph_disconnect = true;
		vector<string> line(mpParameters->mSequenceDegree);
		for (unsigned i = 0; i < mpParameters->mSequenceDegree; ++i) {
			getline(in, line[i]);
		}
		unsigned sequence_length = line[0].length();
		unsigned vertex_counter = 0;
		for (unsigned j = 0; j < sequence_length; j++) {
			for (unsigned i = 0; i < mpParameters->mSequenceDegree; ++i) {
				if (line[i] == "") return;
				char char_label = line[i].at(j);
				string label = stream_cast<string>(char_label);
				if (label == "|") {
					graph_disconnect = true;
					vertex_component_list.push_back(vertex_component);
					vertex_component.clear();
				} else {
					vertex_component.push_back(vertex_counter);
					AddVertexAndEdgesForSequence(oG, label, vertex_counter, graph_disconnect);
					vertex_counter++;
					graph_disconnect = false;
				}
			}
		}
		if (mpParameters->mGraphType == "DIRECTED") AddReverseGraph(oG);
		vertex_component_list.push_back(vertex_component);
		if (vertex_component_list.size() > 1) AddAbstractConnections(oG, vertex_component_list);
	}

	void SetGraphFromSequenceMultiLineTokenFile(istream& in, GraphClass& oG) {
		vector<vector<unsigned> > vertex_component_list;
		vector<unsigned> vertex_component;
		bool graph_disconnect = true;
		vector<vector<string> > multi_line(mpParameters->mSequenceDegree);
		for (unsigned i = 0; i < mpParameters->mSequenceDegree; ++i) {
			vector<string> single_line;
			string line;
			getline(in, line);
			stringstream ss;
			ss << line << endl;
			while (!ss.eof() && ss.good()) {
				string label;
				ss >> label;
				if (label != "") single_line.push_back(label);
			}
			multi_line.push_back(single_line);
		}

		unsigned sequence_length = multi_line[0].size();
		unsigned vertex_counter = 0;
		for (unsigned j = 0; j < sequence_length; j++) {
			for (unsigned i = 0; i < mpParameters->mSequenceDegree; ++i) {
				string label = multi_line[i][j];
				if (label == "|") {
					graph_disconnect = true;
					vertex_component_list.push_back(vertex_component);
					vertex_component.clear();
				} else {
					vertex_component_list.push_back(vertex_component);
					AddVertexAndEdgesForSequence(oG, label, vertex_counter, graph_disconnect);
					vertex_counter++;
					graph_disconnect = false;
				}
			}
		}
		if (mpParameters->mGraphType == "DIRECTED") AddReverseGraph(oG);
		vertex_component_list.push_back(vertex_component);
		if (vertex_component_list.size() > 1) AddAbstractConnections(oG, vertex_component_list);
	}

	void AddVertexAndEdgesForSequence(GraphClass& oG, string aLabel, unsigned aVertexCounter, bool aGraphDisconnect) {
		//set (once) boolean status vectors
		static vector<bool> vertex_status(5, false);
		vertex_status[0] = true; //kernel point
		vertex_status[1] = true; //kind
		vertex_status[2] = true; //viewpoint
		//vertex_status[3] = false; //dead
		//vertex_status[4] = false; //abstraction

		static vector<bool> edge_status(3, false);
		//edge_status[0] = false; //edge dead
		//edge_status[1] = false; //edge abstraction_of
		//edge_status[2] = false; //edge part_of

		unsigned real_vertex_index = oG.InsertVertex();
		vector<string> vertex_symbolic_attribute_list;
		vertex_symbolic_attribute_list.push_back(aLabel);
		oG.SetVertexSymbolicAttributeList(real_vertex_index, vertex_symbolic_attribute_list);
		oG.SetVertexStatusAttributeList(real_vertex_index, vertex_status);
		if (aVertexCounter % mpParameters->mSequenceDegree != 0) {
			oG.SetVertexViewPoint(real_vertex_index, false);
		}
		if (aGraphDisconnect == false) {
			//add edge
			unsigned real_src_index;
			unsigned real_dest_index;
			if (aVertexCounter % mpParameters->mSequenceDegree != 0) { //add edge to preceding vertex
				real_src_index = aVertexCounter;
				real_dest_index = aVertexCounter - 1;
			} else { //vertices with position multiple than mSequenceDegree are connected in sequence
				real_src_index = aVertexCounter;
				real_dest_index = aVertexCounter - mpParameters->mSequenceDegree;
			}
			vector<string> edge_symbolic_attribute_list;
			edge_symbolic_attribute_list.push_back("-"); //NOTE: default edge label is '-'
			unsigned edge_index = oG.InsertEdge(real_src_index, real_dest_index);
			oG.SetEdgeSymbolicAttributeList(edge_index, edge_symbolic_attribute_list);
			oG.SetEdgeStatusAttributeList(edge_index, edge_status);
			if (mpParameters->mGraphType == "UNDIRECTED" || aVertexCounter % mpParameters->mSequenceDegree != 0) { //add reverse edge in case the graph is undirected or for the attribute vertices, so that in both directions these are accessible
				unsigned reverse_edge_index = oG.InsertEdge(real_dest_index, real_src_index);
				oG.SetEdgeSymbolicAttributeList(reverse_edge_index, oG.GetEdgeSymbolicAttributeList(edge_index));
				oG.SetEdgeStatusAttributeList(reverse_edge_index, oG.GetEdgeStatusAttributeList(edge_index));
			}
		}
	}

	void AddAbstractConnections(GraphClass& oG, vector<vector<unsigned> >& aVertexComponentList) {
		//set (once) boolean status vectors
		static vector<bool> vertex_status(5, false);
		//vertex_status[0] = false; //kernel point
		//vertex_status[1] = false; //kind
		//vertex_status[2] = false; //viewpoint
		//vertex_status[3] = false; //dead
		vertex_status[4] = true; //abstraction

		static vector<bool> edge_status(3, false);
		//edge_status[0] = false; //edge dead
		//edge_status[1] = false; //edge abstraction_of
		//edge_status[2] = false; //edge part_of

		//for all pairs of components
		for (unsigned i = 0; i < aVertexComponentList.size(); ++i) {
			for (unsigned j = 0; j < aVertexComponentList.size(); ++j) {
				if (i != j) {
					//join all vertices in one component to all other vertices in the other component
					//add 1 abstract vertex
					unsigned real_vertex_index = oG.InsertVertex();
					vector<string> vertex_symbolic_attribute_list;
					vertex_symbolic_attribute_list.push_back("^L");
					oG.SetVertexSymbolicAttributeList(real_vertex_index, vertex_symbolic_attribute_list);
					oG.SetVertexStatusAttributeList(real_vertex_index, vertex_status);

					for (unsigned ii = 0; ii < aVertexComponentList[i].size(); ii++) { //add part_of edges
						unsigned real_src_index = real_vertex_index;
						unsigned real_dest_index = aVertexComponentList[i][ii];
						vector<string> edge_symbolic_attribute_list;
						edge_symbolic_attribute_list.push_back("@-");
						unsigned edge_index = oG.InsertEdge(real_src_index, real_dest_index);
						oG.SetEdgeSymbolicAttributeList(edge_index, edge_symbolic_attribute_list);
						oG.SetEdgeStatusAttributeList(edge_index, edge_status);
						oG.SetEdgePartOf(edge_index, true);
					}

					for (unsigned jj = 0; jj < aVertexComponentList[j].size(); jj++) { //add abstraction_of edges
						unsigned real_src_index = real_vertex_index;
						unsigned real_dest_index = aVertexComponentList[j][jj];
						vector<string> edge_symbolic_attribute_list;
						edge_symbolic_attribute_list.push_back("^-");
						unsigned edge_index = oG.InsertEdge(real_src_index, real_dest_index);
						oG.SetEdgeSymbolicAttributeList(edge_index, edge_symbolic_attribute_list);
						oG.SetEdgeStatusAttributeList(edge_index, edge_status);
						oG.SetEdgeAbstractionOf(edge_index, true);
					}
				}
			}
		}
	}

	void SetGraphFromGraphGspanFile(istream& in, GraphClass& oG) {
		//status
		vector<bool> vertex_status(5, false);
		vertex_status[0] = true; //kernel point
		vertex_status[1] = true; //kind
		vertex_status[2] = true; //viewpoint
		vertex_status[3] = false; //dead
		vertex_status[4] = false; //abstraction

		vector<bool> edge_status(3, false);
		edge_status[0] = false; //edge dead
		edge_status[1] = false; //edge abstraction_of
		edge_status[2] = false; //edge part_of

		map<string, int> index_map_nominal_to_real;
		string line;
		unsigned line_counter = 0;
		do {
			line_counter++;
			getline(in, line);
			if (line == "") break;
			stringstream ss;
			ss << line << endl;
			char code;
			ss >> code;
			if (code == 't' || code == 'g') break;
			else if (code == 'v' || code == 'V' || code == 'W') {
				//extract vertex id and make map nominal_id -> real_id
				string nominal_vertex_index;
				ss >> nominal_vertex_index;
				unsigned real_vertex_index = oG.InsertVertex();
				index_map_nominal_to_real[nominal_vertex_index] = real_vertex_index;
				//label
				vector<string> vertex_symbolic_attribute_list;
				string label;
				ss >> label;
				vertex_symbolic_attribute_list.push_back(label);
				oG.SetVertexSymbolicAttributeList(real_vertex_index, vertex_symbolic_attribute_list);
				oG.SetVertexStatusAttributeList(real_vertex_index, vertex_status);
				if (code == 'V') oG.SetVertexViewPoint(real_vertex_index, false);
				if (code == 'W') {
					oG.SetVertexViewPoint(real_vertex_index, false);
					oG.SetVertexKernelPoint(real_vertex_index, false);
				}
				char vertex_abstraction_code = label.at(0);
				if (vertex_abstraction_code == '^') {
					oG.SetVertexAbstraction(real_vertex_index, true);
					oG.SetVertexViewPoint(real_vertex_index, false);
					oG.SetVertexKernelPoint(real_vertex_index, false);
				}
			} else if (code == 'e') {
				//extract src and dest vertex id
				string nominal_src_index, nominal_dest_index;
				string label;
				ss >> nominal_src_index >> nominal_dest_index >> label;
				assert(index_map_nominal_to_real.count(nominal_src_index)>0);
				assert(index_map_nominal_to_real.count(nominal_dest_index)>0);
				vector<string> edge_symbolic_attribute_list;
				edge_symbolic_attribute_list.push_back(label);
				unsigned real_src_index = index_map_nominal_to_real[nominal_src_index];
				unsigned real_dest_index = index_map_nominal_to_real[nominal_dest_index];
				unsigned edge_index = oG.InsertEdge(real_src_index, real_dest_index);
				oG.SetEdgeSymbolicAttributeList(edge_index, edge_symbolic_attribute_list);
				oG.SetEdgeStatusAttributeList(edge_index, edge_status);

				char edge_abstraction_code = label.at(0);
				if (edge_abstraction_code == '^') oG.SetEdgeAbstractionOf(edge_index, true);
				if (edge_abstraction_code == '@') oG.SetEdgePartOf(edge_index, true);

				if (mpParameters->mGraphType == "UNDIRECTED" || edge_abstraction_code == '^' || edge_abstraction_code == '@') { //NOTE: edges that are part of the abstraction mechanism should be treated as undirected
					unsigned reverse_edge_index = oG.InsertEdge(real_dest_index, real_src_index);
					oG.SetEdgeSymbolicAttributeList(reverse_edge_index, oG.GetEdgeSymbolicAttributeList(edge_index));
					oG.SetEdgeStatusAttributeList(reverse_edge_index, oG.GetEdgeStatusAttributeList(edge_index));
				}
			} else {
			} //NOTE: ignore other markers
		} while (!in.eof() && in.good());
		if (mpParameters->mGraphType == "DIRECTED") AddReverseGraph(oG);
	}

	void AddReverseGraph(GraphClass& oG) {
		unsigned vsize = oG.VertexSize();
		//add a copy of all vertices
		for (unsigned i = 0; i < vsize; i++) {
			unsigned real_vertex_index = oG.InsertVertex();
			assert(real_vertex_index==i+vsize);
			vector<string> vertex_symbolic_attribute_list = oG.GetVertexSymbolicAttributeList(i);
			for (unsigned t = 0; t < vertex_symbolic_attribute_list.size(); t++) //prepend a prefix to mark the reverse direction
				vertex_symbolic_attribute_list[t] = "r." + vertex_symbolic_attribute_list[t];
			oG.SetVertexSymbolicAttributeList(real_vertex_index, vertex_symbolic_attribute_list);
			oG.SetVertexStatusAttributeList(real_vertex_index, oG.GetVertexStatusAttributeList(i)); //assign original status vector
		}
		//copy all edges swapping src with dest
		for (unsigned i = 0; i < vsize; i++) {
			//get all edges
			vector<unsigned> adj = oG.GetVertexAdjacentList(i);
			for (unsigned j = 0; j < adj.size(); j++) {
				unsigned orig_src = i;
				unsigned orig_dest = adj[j];
				unsigned reverse_src = orig_dest + vsize;
				unsigned reverse_dest = orig_src + vsize;
				unsigned edge_index = oG.InsertEdge(reverse_src, reverse_dest);
				oG.SetEdgeSymbolicAttributeList(edge_index, oG.GetEdgeSymbolicAttributeList(orig_src, orig_dest));
				oG.SetEdgeStatusAttributeList(edge_index, oG.GetEdgeStatusAttributeList(orig_src, orig_dest));
			}
		}
	}

	unsigned Size() {
		return mVectorList.size();
	}
}
;

//------------------------------------------------------------------------------------------------------------------------

class LocalMultiDimensionalScaling {
public:
	Parameters* mpParameters;
	Data* mpData;
	double mTau;
	vector<set<unsigned> > mNeighborhoodList;
	vector<set<unsigned> > mNonNeighborhoodList;
public:
	void Init(Parameters* apParameters, Data* apData) {
		mpParameters = apParameters;
		mpData = apData;
	}

	double Norm(const FVector& aX) {
		return sqrt(dot(aX, aX));
	}

	double Distance(const FVector& aX, const FVector& aZ) {
		FVector diff = combine(aX, 1, aZ, -1);
		return Norm(diff);
	}

	FVector Versor(const FVector& aX, const FVector& aZ) {
		FVector diff = combine(aX, 1, aZ, -1);
		diff.scale(Norm(diff));
		return diff;
	}

	void ComputeLocalMultiDimensionalScaling(vector<FVector>& oXList) {
		double EPSILON = 1 / ((double) mpData->Size() * (double) mpParameters->mLMDSNeighborhoodSize); //NOTE: this corresponds to a change of 10% of the neighbourhood in 5% of the data
		const unsigned NUM_CONSECUTIVE_ITERATIONS_BELOW_EPSILON = 2;
		const unsigned NUM_STEPS_IN_NEIGHBORHOOD_RANGE = 3;
		const double STEP_SIZE_POWER = 1;
		//for random_restart
		vector<FVector> best_x_list;
		double best_distortion = 1;
		double distortion = 1;
		double prev_distortion = 1;
		unsigned best_neighborhood_size = 0;
		double best_log_counter = 0;
		cout << "Computing low dimensional layout normalized distortion" << endl;
		unsigned step = (2 * mpParameters->mLMDSNeighborhoodSizeRange / NUM_STEPS_IN_NEIGHBORHOOD_RANGE);
		step = step < 1 ? 1 : step;
		for (unsigned neighborhood_size_modifier = 0; neighborhood_size_modifier <= 2 * mpParameters->mLMDSNeighborhoodSizeRange; neighborhood_size_modifier += step) {
			unsigned effective_neighborhood_size = mpParameters->mLMDSNeighborhoodSize + neighborhood_size_modifier - mpParameters->mLMDSNeighborhoodSizeRange;
			if (effective_neighborhood_size >= 3) {
				InitNeighbourhoodList(effective_neighborhood_size, mpParameters->mLMDSNonNeighborhoodSize);
				for (double log_counter = 0; log_counter <= mpParameters->mLMDSTauExponentRange; log_counter += STEP_SIZE_POWER) {
					double repulsive_force_tau = mTau * mpParameters->mLMDSTau * pow(10, log_counter);
					for (unsigned random_restart_counter = 0; random_restart_counter < mpParameters->mLMDSNumRandomRestarts; random_restart_counter++) {
						cout
								<< "Neighborhood size: "
								<< effective_neighborhood_size
								<< " ["
								<< mpParameters->mLMDSNeighborhoodSize - mpParameters->mLMDSNeighborhoodSizeRange
								<< ".."
								<< mpParameters->mLMDSNeighborhoodSize + mpParameters->mLMDSNeighborhoodSizeRange
								<< "] "
								<< endl
								<< "Repulsive force: "
								<< repulsive_force_tau
								<< endl
								<< "Restart num "
								<< random_restart_counter + 1
								<< "/"
								<< mpParameters->mLMDSNumRandomRestarts
								<< " (max num iterations: "
								<< mpParameters->mLMDSNumIterations
								<< " or distortion delta less than "
								<< EPSILON
								<< " for more than "
								<< NUM_CONSECUTIVE_ITERATIONS_BELOW_EPSILON
								<< " consecutive iterations)"
								<< endl;

						//initialization phase: random coordinates, computation of long-short distance tradeoff
						vector<FVector> current_x_list;
						LowDimensionalCoordinateInitialization(current_x_list);
						//iterative minimization loop
						ProgressBar pb;
						unsigned low_variation_counter = 0;
						bool low_variation_flag = false;
						for (unsigned itera = 0; itera < mpParameters->mLMDSNumIterations && low_variation_flag == false; ++itera) {
							pb.Count();
							for (unsigned i = 0; i < mpData->Size(); ++i) {
								//forall instances in the neighborhood
								for (set<unsigned>::iterator jt = mNeighborhoodList[i].begin(); jt != mNeighborhoodList[i].end(); ++jt) {
									unsigned j = *jt;
									FVector diff_versor = Versor(current_x_list[j], current_x_list[i]);
									double current_distance = Distance(current_x_list[j], current_x_list[i]);
									double desired_distance = (1 - mpData->ComputeKernel(i, j));
									//double stress = sqrt((current_distance - desired_distance) * (current_distance - desired_distance));
									double stress = fabs(current_distance - desired_distance);
									current_x_list[i].add(diff_versor, stress * mpParameters->mLMDSIterationEpsilon);
								}
								//forall instances not in the neighborhood
								for (set<unsigned>::iterator jt = mNonNeighborhoodList[i].begin(); jt != mNonNeighborhoodList[i].end(); ++jt) {
									unsigned j = *jt;
									FVector diff_versor = Versor(current_x_list[i], current_x_list[j]);
									//double long_distance_repulsive_force = mTau * mpParameters->mLMDSTau * (1 - (double) (itera + 1) / (double) mpParameters->mLMDSNumIterations);
									double long_distance_repulsive_force = repulsive_force_tau;
									current_x_list[i].add(diff_versor, long_distance_repulsive_force * mpParameters->mLMDSIterationEpsilon);
								}
							}
							if (itera % 100 == 0) {
								distortion = Distortion(current_x_list);
								if (fabs(prev_distortion - distortion) < EPSILON) low_variation_counter++;
								else low_variation_counter = 0;
								if (low_variation_counter >= NUM_CONSECUTIVE_ITERATIONS_BELOW_EPSILON) low_variation_flag = true;
								prev_distortion = distortion;
							}
							if (itera % 1000 == 0) {
								ScrambleHighStressInstances(current_x_list);
							}
						}
						distortion = Distortion(current_x_list);
						//keep solution with lowest distortion
						if (distortion < best_distortion) {
							best_x_list = current_x_list;
							best_distortion = distortion;
							best_neighborhood_size = effective_neighborhood_size;
							best_log_counter = log_counter;
							cout << endl << "Saving solution: achieved new low distortion: " << best_distortion << endl;
							SaveEmbedding(best_x_list);
							SaveDistortion(best_x_list);
							SaveNeighbourhoodList();
						} else {
							cout << endl << "Current distortion: " << distortion << endl;
						}
					}
				}
			}
		}
		cout << endl;
		cout << "Best solution found at distortion level: " << best_distortion << " neighbourhood size: " << best_neighborhood_size << " repulsive force multiplicative factor: " << best_log_counter << endl;
		oXList = best_x_list;
	}

	void ComputeLocalMultiDimensionalScalingBatch(vector<FVector>& oXList) {
		double EPSILON = 1 / ((double) mpData->Size() * (double) mpParameters->mLMDSNeighborhoodSize); //NOTE: this corresponds to a change of 10% of the neighbourhood in 5% of the data
		const unsigned NUM_CONSECUTIVE_ITERATIONS_BELOW_EPSILON = 2;
		const unsigned NUM_STEPS_IN_NEIGHBORHOOD_RANGE = 3;
		const double STEP_SIZE_POWER = 1;
		//for random_restart
		vector<FVector> best_x_list;
		double best_distortion = 1;
		double best_stress = 1;
		double distortion = 1;
		double prev_distortion = 1;
		double stress = 1;
		unsigned best_neighborhood_size = 0;
		double best_log_counter = 0;
		cout << "Computing low dimensional layout normalized distortion" << endl;
		unsigned step = (2 * mpParameters->mLMDSNeighborhoodSizeRange / NUM_STEPS_IN_NEIGHBORHOOD_RANGE);
		step = step < 1 ? 1 : step;
		for (unsigned neighborhood_size_modifier = 0; neighborhood_size_modifier <= 2 * mpParameters->mLMDSNeighborhoodSizeRange; neighborhood_size_modifier += step) {
			unsigned effective_neighborhood_size = mpParameters->mLMDSNeighborhoodSize + neighborhood_size_modifier - mpParameters->mLMDSNeighborhoodSizeRange;
			if (effective_neighborhood_size >= 3) {
				InitNeighbourhoodList(effective_neighborhood_size, mpParameters->mLMDSNonNeighborhoodSize);
				for (double log_counter = 0; log_counter <= mpParameters->mLMDSTauExponentRange; log_counter += STEP_SIZE_POWER) {
					double repulsive_force_tau = mTau * mpParameters->mLMDSTau * pow(10, log_counter);
					for (unsigned random_restart_counter = 0; random_restart_counter < mpParameters->mLMDSNumRandomRestarts; random_restart_counter++) {
						cout
								<< "Neighborhood size: "
								<< effective_neighborhood_size
								<< " ["
								<< mpParameters->mLMDSNeighborhoodSize - mpParameters->mLMDSNeighborhoodSizeRange
								<< ".."
								<< mpParameters->mLMDSNeighborhoodSize + mpParameters->mLMDSNeighborhoodSizeRange
								<< "] "
								<< endl
								<< "Repulsive force: "
								<< repulsive_force_tau
								<< endl
								<< "Restart num "
								<< random_restart_counter + 1
								<< "/"
								<< mpParameters->mLMDSNumRandomRestarts
								<< " (max num iterations: "
								<< mpParameters->mLMDSNumIterations
								<< " or distortion delta less than "
								<< EPSILON
								<< " for more than "
								<< NUM_CONSECUTIVE_ITERATIONS_BELOW_EPSILON
								<< " consecutive iterations)"
								<< endl;

						//initialization phase: random coordinates, computation of long-short distance tradeoff
						vector<FVector> current_x_list;
						LowDimensionalCoordinateInitialization(current_x_list);
						//iterative minimization loop
						ProgressBar pb;
						unsigned low_variation_counter = 0;
						bool low_variation_flag = false;
						for (unsigned itera = 0; itera < mpParameters->mLMDSNumIterations && low_variation_flag == false; ++itera) {
							pb.Count();
							vector<FVector> next_x_list = current_x_list;
							for (unsigned i = 0; i < mpData->Size(); ++i) {
								FVector x_delta(mpParameters->mLMDSDimensionality);
								//forall instances in the neighborhood
								for (set<unsigned>::iterator jt = mNeighborhoodList[i].begin(); jt != mNeighborhoodList[i].end(); ++jt) {
									unsigned j = *jt;
									FVector diff_versor = Versor(current_x_list[j], current_x_list[i]);
									double current_distance = Distance(current_x_list[j], current_x_list[i]);
									double desired_distance = (1 - mpData->ComputeKernel(i, j));
									//double stress = sqrt((current_distance - desired_distance) * (current_distance - desired_distance));
									double stress = fabs(current_distance - desired_distance);
									x_delta.add(diff_versor, stress);
								}
								//forall instances not in the neighborhood
								for (set<unsigned>::iterator jt = mNonNeighborhoodList[i].begin(); jt != mNonNeighborhoodList[i].end(); ++jt) {
									unsigned j = *jt;
									FVector diff_versor = Versor(current_x_list[i], current_x_list[j]);
									//double long_distance_repulsive_force = mTau * mpParameters->mLMDSTau * (1 - (double) (itera + 1) / (double) mpParameters->mLMDSNumIterations);
									double long_distance_repulsive_force = repulsive_force_tau;
									x_delta.add(diff_versor, long_distance_repulsive_force);
								}
								x_delta.scale(mpParameters->mLMDSIterationEpsilon);
								next_x_list[i].add(x_delta);
							}
							current_x_list = next_x_list;
							if (itera % 100 == 0) {
								distortion = Distortion(current_x_list);
								if (fabs(prev_distortion - distortion) < EPSILON) low_variation_counter++;
								else low_variation_counter = 0;
								if (low_variation_counter >= NUM_CONSECUTIVE_ITERATIONS_BELOW_EPSILON) low_variation_flag = true;
								prev_distortion = distortion;
							}
							if (itera % 1000 == 0) {
								ScrambleHighStressInstances(current_x_list);
							}
						}
						distortion = Distortion(current_x_list);
						stress = Stress(current_x_list);
						//keep solution with lowest distortion
						if (distortion < best_distortion) {
							best_x_list = current_x_list;
							best_distortion = distortion;
							best_stress = stress;
							best_neighborhood_size = effective_neighborhood_size;
							best_log_counter = log_counter;
							cout << endl << "Saving solution: achieved new low distortion: " << best_distortion << " (stress: " << best_stress << ")" << endl;
							SaveEmbedding(best_x_list);
							SaveDistortion(best_x_list);
							SaveNeighbourhoodList();
						} else {
							cout << endl << "Current distortion: " << distortion << " (stress: " << stress << ")" << endl;
						}
					}
				}
			}
		}
		cout << endl;
		cout << "Best solution found at distortion level: " << best_distortion << " (stress level: " << best_stress << ")" << " neighbourhood size: " << best_neighborhood_size << " repulsive force multiplicative factor: " << best_log_counter << endl;
		oXList = best_x_list;
	}

	void ScrambleHighStressInstances(vector<FVector>& aXList) {
		const double DISTORTION_THRESHOLD = .9;
		const double FRACTION = .9;
		vector<set<unsigned> > low_dim_neighbourhood_list;
		MakeNeighbourhoodList(aXList, low_dim_neighbourhood_list);

		vector<pair<double, unsigned> > stress_list;
		for (unsigned i = 0; i < aXList.size(); ++i) {
			//compute normalized size of the neighbourhood intersection
			set<unsigned> intersection;
			set_intersection(mNeighborhoodList[i].begin(), mNeighborhoodList[i].end(), low_dim_neighbourhood_list[i].begin(), low_dim_neighbourhood_list[i].end(), inserter(intersection, intersection.begin()));
			double distortion = 1 - (double) intersection.size() / sqrt((double) mNeighborhoodList[i].size() * (double) low_dim_neighbourhood_list[i].size());
			if (distortion > DISTORTION_THRESHOLD) {
				//compute stress
				double stress = 0;
				for (unsigned j = 0; j < aXList.size(); ++j) {
					if (i != j) {
						double desired_distance = 1 - mpData->ComputeKernel(i, j);
						double current_distance = Distance(aXList[i], aXList[j]);
						stress += fabs(current_distance - desired_distance);
					}
				}
				double average_stress = stress / (double) (aXList.size() - 1);
				stress_list.push_back(make_pair(-average_stress, i));
			}
		}
		sort(stress_list.begin(), stress_list.end());
		unsigned effective_size = (double) (stress_list.size()) * FRACTION;
		cout << "[" << effective_size << "]";
		for (unsigned i = 0; i < effective_size; ++i) {
			unsigned id = stress_list[i].second;
			FVector x;
			for (unsigned j = 0; j < mpParameters->mLMDSDimensionality; ++j) {
				double value = 2 * random01() - 1;
				x.set(j, value);
			}
			aXList[id] = x;
		}
	}

	void SaveEmbedding(vector<FVector>& aXList) {
		string output_filename = mpParameters->mInputDataFileName + ".embed" + mpParameters->mSuffix;
		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.45: Cannot open file:" + output_filename);
		for (unsigned t = 0; t < aXList.size(); t++) {
			for (unsigned i = 0; i < mpParameters->mLMDSDimensionality; i++) {
				double val = aXList[t].get(i);
				ofs << val << " ";
			}
			ofs << endl;
		}
		cout << "Embedding saved in file " << output_filename << endl;
	}

	void SaveDistortion(vector<FVector>& aXList) {
		string output_filename = mpParameters->mInputDataFileName + ".distortion" + mpParameters->mSuffix;
		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.56: Cannot open file:" + output_filename);

		vector<set<unsigned> > low_dim_neighbourhood_list;
		MakeNeighbourhoodList(aXList, low_dim_neighbourhood_list);
		assert(aXList.size()==mpData->Size());
		for (unsigned i = 0; i < aXList.size(); ++i) {
			//compute normalized size of the neighbourhood intersection
			set<unsigned> intersection;
			set_intersection(mNeighborhoodList[i].begin(), mNeighborhoodList[i].end(), low_dim_neighbourhood_list[i].begin(), low_dim_neighbourhood_list[i].end(), inserter(intersection, intersection.begin()));
			double distortion = 1 - (double) intersection.size() / sqrt((double) mNeighborhoodList[i].size() * (double) low_dim_neighbourhood_list[i].size());
			ofs << distortion << endl;
		}

		cout << "Distortion saved in file " << output_filename << endl;
	}

	void SaveNeighbourhoodList() {
		string output_filename = mpParameters->mInputDataFileName + ".neighbourhood" + mpParameters->mSuffix;
		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.57: Cannot open file:" + output_filename);
		for (unsigned i = 0; i < mpData->Size(); ++i) {
			for (set<unsigned>::iterator it = mNeighborhoodList[i].begin(); it != mNeighborhoodList[i].end(); ++it) {
				ofs << (*it) << " ";
			}
			ofs << endl;
		}
		cout << "Neighbours identity saved in file " << output_filename << endl;
	}

	void InitNeighbourhoodList(unsigned aNeighborhoodSize, unsigned aNonNeighborhoodSize) {
		cout << "Neighbourhood indicators computation" << endl;
		mNeighborhoodList.clear();
		mNonNeighborhoodList.clear();
		//for all instances
		vector<double> all_distance_list;
		{
			ProgressBar pb;
			for (unsigned i = 0; i < mpData->Size(); ++i) {
				vector<pair<double, unsigned> > similarity_list(mpData->Size());
				//determine all pairwise similarities
				for (unsigned j = 0; j < mpData->Size(); ++j) {
					double distance = 1 - mpData->ComputeKernel(i, j);
					similarity_list[j] = make_pair(distance, j);
					all_distance_list.push_back(distance);
				}
				//sort and retain the closest k indices
				sort(similarity_list.begin(), similarity_list.end());
				set<unsigned> neighbour_list;
				set<unsigned> non_neighbour_list;
				unsigned effective_size = min(aNeighborhoodSize, (unsigned) similarity_list.size());
				for (unsigned k = 0; k < effective_size; ++k) {
					unsigned neighbour_id = similarity_list[k].second;
					if (neighbour_id != i) neighbour_list.insert(neighbour_id);
				}
				unsigned step = (similarity_list.size() - 2 * effective_size) / aNonNeighborhoodSize;
				step = (step == 0) ? 1 : step;
				for (unsigned k = 2 * effective_size; k < similarity_list.size(); k = k + step) {
					unsigned neighbour_id = similarity_list[k].second;
					non_neighbour_list.insert(neighbour_id);
				}
				mNeighborhoodList.push_back(neighbour_list);
				mNonNeighborhoodList.push_back(non_neighbour_list);
				pb.Count();
			}
		}
		//init Tau
		cout << "Computing local-global tradeoff factor: ";
		sort(all_distance_list.begin(), all_distance_list.end());
		unsigned median_index = all_distance_list.size() / 2;
		double median_distance = all_distance_list[median_index];
		mTau = (double) mNeighborhoodList[0].size() / (double) mNonNeighborhoodList[0].size() * median_distance;
		cout << mTau << endl;
	}

	void MakeNeighbourhoodList(vector<FVector>& aXList, vector<set<unsigned> >& oNeighbourhoodList) {
		for (unsigned i = 0; i < aXList.size(); ++i) {
			vector<pair<double, unsigned> > similarity_list(aXList.size());
			//determine all pairwise similarities
			for (unsigned j = 0; j < aXList.size(); ++j) {
				if (i != j) { //NOTE: exclude self
					FVector difference_x = combine(aXList[i], 1, aXList[j], -1);
					double distance = dot(difference_x, difference_x);
					similarity_list[j] = make_pair(distance, j);
				}
			}
			//sort and retain the closest k indices
			sort(similarity_list.begin(), similarity_list.end());
			set<unsigned> neighbour_list;
			unsigned effective_size = min(mpParameters->mLMDSNeighborhoodSize, (unsigned) similarity_list.size());
			for (unsigned k = 0; k < effective_size; ++k) {
				unsigned neighbour_id = similarity_list[k].second;
				neighbour_list.insert(neighbour_id);
			}
			oNeighbourhoodList.push_back(neighbour_list);
		}
	}

	void LowDimensionalCoordinateInitialization(vector<FVector>& oXList) {
		unsigned choice = randomUnsigned(2);
		choice = 1;
		if (choice == 0) KernelLowDimensionalCoordinateInitialization(oXList);
		else RandomLowDimensionalCoordinateInitialization(oXList);
	}

	void KernelLowDimensionalCoordinateInitialization(vector<FVector>& oXList) {
		oXList.clear();
		//select mLMDSDimensionality random instances
		vector<unsigned> anchor_list;
		//Randomly permute consecutive index sequence
		for (unsigned i = 0; i < mpData->Size(); ++i)
			anchor_list.push_back(i);
		PermuteVector<unsigned>(anchor_list);
		anchor_list.resize(mpParameters->mLMDSDimensionality);
		cout << "Kernel based initialization. ";
		cout << "Selected reference instances: ";
		for (unsigned i = 0; i < anchor_list.size(); ++i)
			cout << anchor_list[i] << " ";
		cout << endl;

		//set the coordinates as the distance from those instances
		for (unsigned t = 0; t < mpData->Size(); ++t) {
			FVector x;
			for (unsigned k = 0; k < mpParameters->mLMDSDimensionality; ++k) {
				double distance = 1 - mpData->ComputeKernel(anchor_list[k], t);
				x.set(k, distance);
			}
			oXList.push_back(x);
		}
	}

	void RandomLowDimensionalCoordinateInitialization(vector<FVector>& oXList) {
		cout << "Random initialization. " << endl;
		oXList.clear();
		const double span = 2;
		for (unsigned i = 0; i < mpData->Size(); ++i) {
			FVector x;
			for (unsigned j = 0; j < mpParameters->mLMDSDimensionality; ++j) {
				double value = span * random01() - span / 2;
				x.set(j, value);
			}
			oXList.push_back(x);
		}
	}

	double Distortion(vector<FVector>& aXList) {
		unsigned local_continuity = 0;
		vector<set<unsigned> > low_dim_neighbourhood_list;
		MakeNeighbourhoodList(aXList, low_dim_neighbourhood_list);
		assert(aXList.size()==mpData->Size());
		for (unsigned i = 0; i < aXList.size(); ++i) {
			//compute size of the neighborhood intersection
			set<unsigned> intersection;
			set_intersection(mNeighborhoodList[i].begin(), mNeighborhoodList[i].end(), low_dim_neighbourhood_list[i].begin(), low_dim_neighbourhood_list[i].end(), inserter(intersection, intersection.begin()));
			local_continuity += intersection.size();
		}
		double average_local_continuity = ((double) local_continuity / (double) mpParameters->mLMDSNeighborhoodSize) / (double) mpData->Size();
		double average_local_continuity_adjusted_for_chance = average_local_continuity - (double) mpParameters->mLMDSNeighborhoodSize / (double) mpData->Size();
		return 1 - average_local_continuity_adjusted_for_chance;
	}

	double Stress(vector<FVector>& aXList) {
		assert(aXList.size()==mpData->Size());
		double stress = 0;

		//compute average difference of distances
		unsigned counter = 0;
		for (unsigned i = 0; i < aXList.size(); ++i) {
			for (unsigned j = 0; j < aXList.size(); ++j) {
				if (i != j) {
					double desired_distance = 1 - mpData->ComputeKernel(i, j);
					double current_distance = Distance(aXList[i], aXList[j]);
					stress += fabs(current_distance - desired_distance) / (desired_distance);
					//stress += sqrt((current_distance - desired_distance) * (current_distance - desired_distance))/(desired_distance) ;
					counter++;
				}
			}
		}
		double average_stress = stress / (double) counter;
		return average_stress;
	}

};

//------------------------------------------------------------------------------------------------------------------------
/**
 Encapsulates a linear SVM model trainable with stochastic gradient
 descent over graph instances explicitly mapped by the NSPDK kernel
 */
class CoreModel {
	/**
	 Data structure to: 1) facilitate rebalancing of dataset by copying
	 multiple times a reference to the instance; and 2) retrieve
	 prediction efficiently by overwriting a reference to the margin
	 list cell element.
	 */
	struct TrainItem {
		int mInstanceID;
		double mTarget;
		double* mpMargin;
		SVector* mpInstance;
	};
public:
	Parameters* mpParameters;
protected:
	double mLambda;
	unsigned mEpochs;
	double mWScale;
	double mBias;
	SVector mW;
public:
	CoreModel() {
	}

	void Init(Parameters* apParameters) {
		mpParameters = apParameters;
		mLambda = apParameters->mLambda;
		mEpochs = apParameters->mEpochs;
		mWScale = 1;
		mBias = 0;
	}

	/**
	 Prints several informative measures given a list of predictions and a list of true targets
	 */
	double OutputPerformanceMeasures(ostream& out, const vector<double>& aMarginList, const vector<double>& aTargetList) {
		assert(aMarginList.size()==aTargetList.size());
		unsigned size = aMarginList.size();
		unsigned error = 0;
		unsigned correct = 0;
		unsigned tp, tn, fp, fn;
		tp = tn = fp = fn = 0;
		for (unsigned i = 0; i < aMarginList.size(); ++i) {
			double margin = aMarginList[i];
			double prediction = margin > 0 ? 1 : -1;
			double target = aTargetList[i];
			if (prediction != target) error++;
			if (prediction == target) correct++;
			if (prediction > 0 && target > 0) tp++;
			if (prediction > 0 && target < 0) fp++;
			if (prediction < 0 && target > 0) fn++;
			if (prediction < 0 && target < 0) tn++;
		}

		double pprecision = (double) tp / (tp + fp);
		double precall = (double) tp / (tp + fn);
		double pfmeasure = 2 * pprecision * precall / (pprecision + precall);

		double nprecision = (double) tn / (tn + fn);
		double nrecall = (double) tn / (tn + fp);
		double nfmeasure = 2 * nprecision * nrecall / (nprecision + nrecall);

		double bprecision = (pprecision + nprecision) / 2;
		double brecall = (precall + nrecall) / 2;
		double bfmeasure = (pfmeasure + nfmeasure) / 2;

		out << TAB << "Size: " << size << endl;
		out << TAB << "Correct: " << correct << " ( " << correct * 100 / (double) (correct + error) << " %)" << endl;
		out << TAB << "Error: " << error << " ( " << error * 100 / (double) (correct + error) << " %)" << endl;
		out << TAB << "Confusion table:" << endl;
		out << TAB << "TP:" << tp << " FP:" << fp << endl;
		out << TAB << "FN:" << fn << " TN:" << tn << endl;
		out << TAB << "+Precision:" << pprecision << " +Recall:" << precall << " +F-measure:" << pfmeasure << endl;
		out << TAB << "-Precision:" << nprecision << " -Recall:" << nrecall << " -F-measure:" << nfmeasure << endl;
		out << TAB << "bPrecision:" << bprecision << " bRecall:" << brecall << " bF-measure:" << bfmeasure << endl;

		return bfmeasure;
	}

	void Save(ostream& out) {
		out << "bias " << mBias << endl;
		out << "wscale " << mWScale << endl;
		out << "w " << mW;
	}

	void Load(istream& in) {
		string attribute = "";
		string expected = "";
		in >> attribute >> mBias;
		expected = "bias";
		if (attribute != expected) throw range_error("ERROR2.17: Format error: expecting [" + expected + "] but found [" + attribute + "]");
		in >> attribute >> mWScale;
		expected = "wscale";
		if (attribute != expected) throw range_error("ERROR2.18: Format error: expecting [" + expected + "] but found [" + attribute + "]");
		in >> attribute >> mW;
		assert(attribute=="w");
	}

	vector<double> Train(vector<double>& aTargetList, vector<unsigned>& aTrainsetIDList, Data& aData) {
		if (aTrainsetIDList.size() != aTargetList.size()) throw range_error("ERROR2.19: Data list and Target list have not the same size: #data:" + stream_cast<string>(aTrainsetIDList.size()) + " #targets:" + stream_cast<string>(aTargetList.size()));

		vector<SVector*> sv_data_list;
		for (unsigned i = 0; i < aTrainsetIDList.size(); ++i) {
			unsigned id = aTrainsetIDList[i];
			sv_data_list.push_back(&aData.mVectorList[id]);
		}

		//allocate local target list and compute positive/negative target counts
		vector<double> target_list;
		unsigned p, n;
		p = n = 0;
		for (unsigned i = 0; i < aTargetList.size(); ++i) {
			if (aTargetList[i] > 0) p++;
			else n++;
			target_list.push_back(aTargetList[i]);
		}

		//if no instance has negative class then generate negative instances with opposite features wrt positive instances
		vector<SVector> synth_neg_sv_data_list(aTrainsetIDList.size());
		if (n == 0) {
			cout << "No negative instances: proceeding to generate " << aTrainsetIDList.size() << " negative instances with opposite features wrt positive instances" << endl;
			ProgressBar pb;
			for (unsigned i = 0; i < aTrainsetIDList.size(); ++i) {
				SVector x = *sv_data_list[i];
				x.scale(-1);
				synth_neg_sv_data_list[i] = x;
				sv_data_list.push_back(&(synth_neg_sv_data_list[i]));
				target_list.push_back(-1);
				pb.Count();
			}
		}

		//clear margin list and allocate memory
		vector<double> margin_list(target_list.size());
		//...rebalance classes use pointers array to scramble and oversample
		vector<TrainItem> balanced_dataset;
		BalanceDataset(aTrainsetIDList, target_list, margin_list, sv_data_list, balanced_dataset);

		//train on balanced train data
		CoreTrainRoutine(balanced_dataset);

		//output statistics on original train data (no class balance)
		OutputModelInfo();
		cout << "Performance on train set:" << endl;
		OutputPerformanceMeasures(cout, margin_list, target_list);

		if (n == 0) { //if no instance has negative class then re-test model only on real training data to extract margins and predictions
			margin_list = Test(aTrainsetIDList, aData);
		}
		return margin_list;
	}

	void BalanceDataset(vector<unsigned>& aDatasetIDList, vector<double>& aTargetList, vector<double>& oMarginList, vector<SVector*>& aSVDataList, vector<TrainItem>& oDataset) {
		//compute class distribution
		unsigned p, n;
		p = n = 0;
		for (unsigned i = 0; i < aTargetList.size(); ++i)
			if (aTargetList[i] > 0) p++;
			else n++;
		cout << "Class distribution: " << p + n << " (+:" << p << " -:" << n << ") " << "[+:" << (double) p / (p + n) << " -:" << (double) n / (p + n) << "]" << endl;

		//separate positive from negative instances
		vector<TrainItem> positive_data_list;
		vector<TrainItem> negative_data_list;
		if (aTargetList.size() != aSVDataList.size()) throw range_error("ERROR2.20: number of target values: " + stream_cast<string>(aTargetList.size()) + " is different from dataset size:" + stream_cast<string>(aSVDataList.size()));
		for (unsigned i = 0; i < aTargetList.size(); ++i) {
			TrainItem ti;
			ti.mTarget = aTargetList[i];
			ti.mpInstance = aSVDataList[i];
			ti.mpMargin = &oMarginList[i];
			if (i < aDatasetIDList.size()) { //Synthesized instances are appended after the real instances, so the size information of the original id_list marks the start of the syntesized instances
				ti.mInstanceID = aDatasetIDList[i];
			} else ti.mInstanceID = -1; //if the instance has been synthesized then it has no correspondent original graph
			if (aTargetList[i] == 1) positive_data_list.push_back(ti);
			else if (aTargetList[i] == -1) negative_data_list.push_back(ti);
			else throw range_error("ERROR2.21: target has to be 1 or -1; cannot be: " + stream_cast<string>(aTargetList[i]));
		}
		//randomly shuffle data
		for (unsigned i = 0; i < positive_data_list.size(); ++i) {
			unsigned j = rand() * positive_data_list.size() / RAND_MAX;
			swap(positive_data_list[i], positive_data_list[j]);
		}
		for (unsigned i = 0; i < negative_data_list.size(); ++i) {
			unsigned j = rand() * negative_data_list.size() / RAND_MAX;
			swap(negative_data_list[i], negative_data_list[j]);
		}

		//over-sample minority class only if there is an imbalance higher than MIN_KFOLD_IMBALANCE and if there is at least one instance for the minority class
		vector<TrainItem> balanced_positive_data_list;
		vector<TrainItem> balanced_negative_data_list;
		const double MIN_KFOLD_IMBALANCE = 1;
		if (p != 0 && p < n / MIN_KFOLD_IMBALANCE) {
			cout << "Oversampling positive factor: " << n / (double) p << endl;
			unsigned ratio = n / p;
			unsigned reminder = n % p;
			//duplicate a number of times equal to ratio the datastaset itself
			for (unsigned i = 0; i < ratio; i++)
				balanced_positive_data_list.insert(balanced_positive_data_list.end(), positive_data_list.begin(), positive_data_list.end());
			//add the remainder instances
			for (unsigned i = 0; i < reminder; i++)
				balanced_positive_data_list.push_back(positive_data_list[i]);
			balanced_negative_data_list = negative_data_list;
		} else if (n != 0 && n < p / MIN_KFOLD_IMBALANCE) {
			cout << "Oversampling negative factor: " << p / (double) n << endl;
			unsigned ratio = p / n;
			unsigned reminder = p % n;
			for (unsigned i = 0; i < ratio; i++)
				balanced_negative_data_list.insert(balanced_negative_data_list.end(), negative_data_list.begin(), negative_data_list.end());
			for (unsigned i = 0; i < reminder; i++)
				balanced_negative_data_list.push_back(negative_data_list[i]);
			balanced_positive_data_list = positive_data_list;
		} else {
			balanced_positive_data_list = positive_data_list;
			balanced_negative_data_list = negative_data_list;
		}

		//compose dataset by alternating positive and negative examples
		unsigned i;
		for (i = 0; i < balanced_positive_data_list.size(); i++) {
			oDataset.push_back(balanced_positive_data_list[i]);
			if (i < balanced_negative_data_list.size()) oDataset.push_back(balanced_negative_data_list[i]);
		}
		for (unsigned j = i; j < balanced_negative_data_list.size(); j++)
			oDataset.push_back(balanced_negative_data_list[i]);

		//compute new class ratio
		unsigned bp = 0, bn = 0;
		for (unsigned i = 0; i < oDataset.size(); i++)
			if (oDataset[i].mTarget > 0) bp++;
			else bn++;
		cout << "Rebalanced dataset: " << bp + bn << " (+:" << bp << " -:" << bn << ")" << endl;
	}

	void CoreTrainRoutine(vector<TrainItem>& aDataset) {
#ifdef DEBUGON
		VectorClass s;
#endif
		// Shift t in order to have a reasonable initial learning rate. This assumes |x| \approx 1.
		double maxw = 1.0 / sqrt(mLambda);
		double typw = sqrt(maxw);
		double eta0 = typw / max(1.0, dloss(-typw));
		double t = 1 / (eta0 * mLambda);

		//Iterate epochs times in gradient descent
		ProgressBar pb(1);
		OutputTrainingInfo();
		cout << "Training for " << mEpochs << " epochs." << endl;
		for (unsigned e = 0; e < mEpochs; e++) {
			pb.Count();
			//iterate over all train instances
			for (unsigned i = 0; i < aDataset.size(); ++i) {
				double eta = 1.0 / (mLambda * t);
				double s = 1 - eta * mLambda;
				mWScale *= s;
				if (mWScale < 1e-9) {
					mW.scale(mWScale);
					mWScale = 1;
				}
				const SVector &x = (*aDataset[i].mpInstance);
				double y = aDataset[i].mTarget;
				double wx = dot(mW, x) * mWScale;
				double margin = (wx + mBias);
				(*(aDataset[i].mpMargin)) = margin;
				double z = y * margin;
#if LOSS < LOGLOSS
				if (z < 1)
#endif
						{
					double etd = eta * dloss(z);
					mW.add(x, etd * y / mWScale);
#if BIAS
					// Slower rate on the bias because it learns at each iteration.
					mBias += etd * y * 0.01;
#endif
				}
				t += 1;
			}
#ifdef DEBUGON
			s.PushBack(mW.sparse_size());
#endif
		}
#ifdef DEBUGON
		cout << endl << "W size statistics: ";
		s.OutputStatistics(cout);
		cout << endl;
#endif
	}

	vector<double> Test(vector<unsigned>& aTestSetIDList, Data& aData) {
		vector<double> margin_list;
		for (unsigned i = 0; i < aTestSetIDList.size(); ++i) {
			unsigned gid = aTestSetIDList[i];
			SVector& x = aData.mVectorList[gid];
			double y = Predict(x);
			margin_list.push_back(y);
		}
		return margin_list;
	}

	inline double Predict(const SVector& x) {
		return dot(mW, x) * mWScale + mBias;
	}

	void OutputTrainingInfo() {
		cout << SEP << endl;
		cout << "Training information" << endl;
		cout << SEP << endl;
		cout << "Lambda: " << mLambda << endl;
		cout << "Epochs: " << mEpochs<< endl;
		cout << SEP << endl;
	}

	void OutputModelInfo() {
		cout << SEP << endl;
		cout << "Model information" << endl;
		cout << SEP << endl;
		cout << "W Norm: " << sqrt(dot(mW, mW)) * mWScale << endl;
		cout << "Bias: " << mBias << endl;
		cout << SEP << endl;
	}

};

//------------------------------------------------------------------------------------------------------------------------
class Model {
protected:
	Parameters mParameters;
	Data mData;
	CoreModel mCoreModel;

public:
	Model() {
	}
	void Init(int argc, const char **argv) {
		mParameters.Init(argc, argv);
		srand(mParameters.mRandomSeed);
		mData.Init(&mParameters);
		mCoreModel.Init(&mParameters);
	}
	void Exec() {
		switch (mParameters.mActionCode) {
		case TRAIN:
			mData.LoadTarget();
			mData.LoadData();
			Train();
			break;
		case CROSS_VALIDATION: {
			mData.LoadTarget();
			mData.LoadData();
			map<unsigned, vector<double> > results;
			CrossValidation(results);
		}
			break;
		case CONFIDENCE:
			mData.LoadTarget();
			mData.LoadData();
			Confidence();
			break;
		case PARAMETERS_OPTIMIZATION:
			mData.LoadTarget();
			ParametersOptimization();
			break;
		case LEARNING_CURVE:
			mData.LoadTarget();
			mData.LoadData();
			LearningCurve();
			break;
		case MATRIX:
			mData.LoadData();
			Matrix();
			break;
		case EMBED:
			mData.LoadData();
			ComputeLocalMultiDimensionalScaling();
			break;
		case TEST:
		case TEST_PART:
		case FEATURE:
		case FEATURE_PART:
			OnlineProcessing();
			break;
		default:
			throw range_error("ERROR2.2: Unknown action parameter: " + mParameters.mAction);
		}
	}

	void OnlineProcessing() {
		mData.mKernel.ParametersSetup();
		igzstream fin;
		fin.open(mParameters.mInputDataFileName.c_str());
		if (!fin) throw range_error("ERROR2.11: Cannot open file: " + mParameters.mInputDataFileName);

		//compose output filename
		string output_filename = mParameters.mInputDataFileName;
		switch (mParameters.mActionCode) {
		case TEST:
			output_filename += ".prediction";
			break;
		case TEST_PART:
			output_filename += ".prediction_part";
			break;
		case FEATURE:
			output_filename += ".feature";
			break;
		case FEATURE_PART:
			output_filename += ".feature_part";
			break;
		default:
			throw range_error("ERROR2.2: Unknown action parameter: " + mParameters.mAction);
		}
		output_filename += mParameters.mSuffix;

		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.16: Cannot open file: " + output_filename);

		//init phase
		switch (mParameters.mActionCode) {
		case TEST:
		case TEST_PART:
			LoadModel();
			break;
		case FEATURE:
		case FEATURE_PART:
			//no initialization needed
			break;
		default:
			throw range_error("ERROR2.2: Unknown action parameter: " + mParameters.mAction);
		}

		//perform online action
		cout << "Loading file:" << mParameters.mInputDataFileName << endl;
		{
			ProgressBar pb;
			unsigned instance_counter = 0;
			while (!fin.eof()) {
				GraphClass g;
				mData.SetGraphFromFile(fin, g);
				if (!g.IsEmpty()) {
					switch (mParameters.mActionCode) {
					case TEST: {
						SVector x;
						mData.mKernel.GenerateFeatureVector(g, x);
						double margin = mCoreModel.Predict(x);
						int prediction = margin > 0 ? 1 : -1;
						ofs << prediction << " " << margin << endl;
					}
						break;
					case TEST_PART: {
						vector<SVector> graph_vertex_vector_list;
						mData.mKernel.GenerateVertexFeatureVector(g, graph_vertex_vector_list);
						//for each vertex, compute margin
						unsigned size = mParameters.mGraphType == "DIRECTED" ? graph_vertex_vector_list.size() / 2 : graph_vertex_vector_list.size();
						for (unsigned vertex_id = 0; vertex_id < size; ++vertex_id) {
							double margin = mCoreModel.Predict(graph_vertex_vector_list[vertex_id]);
							if (mParameters.mGraphType == "DIRECTED") margin += mCoreModel.Predict(graph_vertex_vector_list[vertex_id + size]);
							ofs << instance_counter << " " << vertex_id << " " << margin << endl;
						}
					}
						break;
					case FEATURE: {
						SVector x;
						mData.mKernel.GenerateFeatureVector(g, x);
						ofs << x;
					}
						break;
					case FEATURE_PART: {
						vector<SVector> graph_vertex_vector_list;
						mData.mKernel.GenerateVertexFeatureVector(g, graph_vertex_vector_list);
						unsigned size = mParameters.mGraphType == "DIRECTED" ? graph_vertex_vector_list.size() / 2 : graph_vertex_vector_list.size();
						for (unsigned vertex_id = 0; vertex_id < size; ++vertex_id) {
							SVector x = graph_vertex_vector_list[vertex_id];
							if (mParameters.mGraphType == "DIRECTED") x.add(graph_vertex_vector_list[vertex_id + size]);
							ofs << instance_counter << " " << vertex_id << " " << x;
						}
					}
						break;
					default:
						throw range_error("ERROR2.2: Unknown action parameter: " + mParameters.mAction);
					}
					pb.Count();
					instance_counter++;
				}
			}
		}
		cout << "Result saved in file " + output_filename << endl;
	}

	void Train() {
		cout << endl << SEP << endl << "Train phase" << endl << SEP << endl;
		ProgressBar pb;
		pb.Count();

		vector<unsigned> train_id_list;
		for (unsigned i = 0; i < mData.mTargetList.size(); ++i)
			train_id_list.push_back(i);
		Train(mData.mTargetList, train_id_list, mData);
		SaveModel();

#ifdef DEBUGON
		string feature_map_file_name = mParameters.mInputDataFileName + "_feature.map" + mParameters.mSuffix;
		cout << endl << "Saving feature map to file: " << feature_map_file_name << endl;
		ofstream ofs(feature_map_file_name.c_str());
		mData.mKernel.mpFeatureGenerator->OutputFeatureMap(ofs);
		ofs.close();
#endif
		cout << endl << "Train phase completed:";
	}

	void ParametersOptimization() {
		string output_filename = mParameters.mInputDataFileName + ".opt_param" + mParameters.mSuffix;
		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.45: Cannot open file:" + output_filename);

		ProgressBar pbt;
		pbt.Count();
		double current_bfmeasure = 0;
		double best_bfmeasure = 0;
		unsigned r_max = mParameters.mRadius;
		unsigned d_max = mParameters.mDistance;
		double lambda_max = mParameters.mLambda;
		unsigned epochs_max = mParameters.mEpochs;

		//default values
		mParameters.mRadius = 1;
		mParameters.mDistance = 1;
		mParameters.mLambda = 1e-4;
		mParameters.mEpochs = 5;
		mData.LoadData();

		const unsigned NUM_LINE_SEARCH_ITERATIONS = 3;
		cerr << "Line search parameters optimization: " << NUM_LINE_SEARCH_ITERATIONS << " iterations." << endl;
		cerr << "Initial parameters configuration: -r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;

		for (unsigned line_search_iteration = 0; line_search_iteration < NUM_LINE_SEARCH_ITERATIONS; line_search_iteration++) {
			cerr << "Iteration " << line_search_iteration + 1 << "/" << NUM_LINE_SEARCH_ITERATIONS << endl;

			//optimize lambda
			best_bfmeasure = 0;
			double lambda_best = mParameters.mLambda;
			const double LAMBDA_UPPER_BOUND = 0.01;
			//double lambda_step = exp((log(LAMBDA_UPPER_BOUND) - log(lambda_max)) / (2 * mParameters.mLearningCurveNumPoints));
			double lambda_step = 10;
			for (double lambda = lambda_max; lambda <= LAMBDA_UPPER_BOUND; lambda *= lambda_step) {
				mParameters.mLambda = lambda;
				map<unsigned, vector<double> > results;
				current_bfmeasure = CrossValidation(results);
				if (current_bfmeasure > best_bfmeasure) {
					best_bfmeasure = current_bfmeasure;
					lambda_best = lambda;
					cerr << "*";
				}
				cerr << "l:" << lambda << " bf:" << current_bfmeasure << " " << endl;
			}
			mParameters.mLambda = lambda_best;
			cerr << "Lambda optimization: bFmeasure:" << best_bfmeasure << " current best parameters configuration: -r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;

			//optimize epochs
			best_bfmeasure = 0;
			unsigned epochs_best = mParameters.mEpochs;
			unsigned epochs_step = (double) epochs_max / (double) mParameters.mLearningCurveNumPoints;
			epochs_step=epochs_step==0?1:epochs_step;
			for (unsigned epochs = epochs_step; epochs <= epochs_max; epochs += epochs_step) {
				mParameters.mEpochs = epochs;
				map<unsigned, vector<double> > results;
				current_bfmeasure = CrossValidation(results);
				if (current_bfmeasure > best_bfmeasure) {
					best_bfmeasure = current_bfmeasure;
					epochs_best = epochs;
					cerr << "*";
				}
				cerr << "e:" << epochs << " bf:" << current_bfmeasure << " " << endl;
			}
			mParameters.mEpochs = epochs_best;
			cerr << "Epochs optimization: bFmeasure:" << best_bfmeasure << " current best parameters configuration: -r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;

			//optimize distance
			best_bfmeasure = 0;
			unsigned d_best = mParameters.mDistance;
			for (unsigned d = 0; d <= d_max; d++) {
				mParameters.mDistance = d;
				mData.LoadData();
				map<unsigned, vector<double> > results;
				current_bfmeasure = CrossValidation(results);
				if (current_bfmeasure > best_bfmeasure) {
					best_bfmeasure = current_bfmeasure;
					d_best = d;
					cerr << "*";
				}
				cerr << "d:" << d << " bf:" << current_bfmeasure << " " << endl;
			}
			mParameters.mDistance = d_best;
			cerr << "Distance optimization: bFmeasure:" << best_bfmeasure << " current best parameters configuration: -r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;

			//optimize radius
			best_bfmeasure = 0;
			unsigned r_best = mParameters.mRadius;
			for (unsigned r = 0; r <= r_max; r++) {
				mParameters.mRadius = r;
				mData.LoadData();
				map<unsigned, vector<double> > results;
				current_bfmeasure = CrossValidation(results);
				if (current_bfmeasure > best_bfmeasure) {
					best_bfmeasure = current_bfmeasure;
					r_best = r;
					cerr << "*";
				}
				cerr << "r:" << r << " bf:" << current_bfmeasure << " " << endl;
			}
			mParameters.mRadius = r_best;
			cerr << "Radius optimization: bFmeasure:" << best_bfmeasure << " current best parameters configuration: -r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;
		}
		cerr << "Best parameter configuration obtains a balanced F-measure= " << best_bfmeasure << endl;
		cerr << "For best predictive performance use the following parameter setting:" << endl;
		cerr << "-r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;
		ofs << "-r " << mParameters.mRadius << " -d " << mParameters.mDistance << " -e " << mParameters.mEpochs << " -l " << mParameters.mLambda << endl;
		cout << "Optimal parameters configuration saved in file " << output_filename << endl;
	}

	void Confidence() {
		string output_filename = mParameters.mInputDataFileName + ".conf" + mParameters.mSuffix;
		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.45: Cannot open file:" + output_filename);

		ProgressBar pbt;
		pbt.Count();
		vector<map<unsigned, vector<double> > > multi_result_list;
		for (unsigned i = 0; i < mParameters.mLearningCurveNumPoints; ++i) {
			mParameters.mRandomSeed += i;
			map<unsigned, vector<double> > result_list;
			CrossValidation(result_list);
			multi_result_list.push_back(result_list);
		}

		unsigned size = mData.Size();
		vector<double> accuracy_list(size);
		vector<VectorClass> margin_list(size);
		for (unsigned r = 0; r < mParameters.mLearningCurveNumPoints; ++r) {
			for (unsigned i = 0; i < size; ++i) {
				double target = multi_result_list[r][i][0];
				double prediction = multi_result_list[r][i][1];
				double margin=multi_result_list[r][i][2];
				margin_list[i].PushBack(margin);
				if (prediction == target) accuracy_list[i] += 1 / (double) mParameters.mLearningCurveNumPoints;
			}
		}

		//write results to file
		for (unsigned i = 0; i < size; ++i) {
			double target=multi_result_list[0][i][0];
			double mean_margin=margin_list[i].Mean();
			double std_margin=margin_list[i].StandardDeviation();
			ofs << target<<" "<<accuracy_list[i] << " "<<mean_margin <<" "<<std_margin<< endl;
		}
		cout<< "Target, average accuracy, average margin, margin standard deviation saved in file " << output_filename << endl;
	}

	double CrossValidation(map<unsigned, vector<double> >& oTestResultMap) {
		oTestResultMap.clear();
		srand(mParameters.mRandomSeed);
		ProgressBar pbt;
		pbt.Count();
		vector<pair<double, double> > prediction_list;
		vector<pair<double, double> > margin_list;

		//randomly shuffle indices
		unsigned size = mData.Size();
		vector<unsigned> data_id_list;
		for (unsigned i = 0; i < size; ++i)
			data_id_list.push_back(i);
		for (unsigned i = 0; i < size; ++i) {
			unsigned j = rand() * size / RAND_MAX;
			swap(data_id_list[i], data_id_list[j]);
		}

		//loop to build train-test split in the cross validation way
		for (unsigned f = 0; f < mParameters.mCrossValidationNumFolds; f++) {
			ProgressBar pbf;
			pbf.Count();
			cout << SEP << endl;
			cout << TAB << TAB << "Fold: " << f + 1 << " of " << mParameters.mCrossValidationNumFolds << endl;
			vector<unsigned> train_id_list;
			vector<unsigned> test_id_list;
			map<unsigned, unsigned> class_counter_map;
			for (unsigned i = 0; i < size; ++i) {
				unsigned id = data_id_list[i];
				double target = mData.mTargetList[id];
				if (class_counter_map.count(target) == 0) class_counter_map[target] = 1;
				else class_counter_map[target]++;
				if (target != 0 && class_counter_map[target] % mParameters.mCrossValidationNumFolds == f) //NOTE: exclude un supervised material from test element list
				test_id_list.push_back(id);
				else train_id_list.push_back(id);
			}
			//sort the indices in order to guarantee sequential file access
			sort(train_id_list.begin(), train_id_list.end());
			sort(test_id_list.begin(), test_id_list.end());
			//extract target list for training
			vector<double> train_target_list;
			for (unsigned i = 0; i < train_id_list.size(); i++) {
				unsigned id = train_id_list[i];
				train_target_list.push_back(mData.mTargetList[id]);
			}
			//extract target list for testing
			vector<double> test_target_list;
			for (unsigned i = 0; i < test_id_list.size(); i++) {
				unsigned id = test_id_list[i];
				test_target_list.push_back(mData.mTargetList[id]);
			}

			//perform training
			Train(train_target_list, train_id_list, mData);
			string model_filename_prefix = "_" + stream_cast<string>(f + 1);
			SaveModel(model_filename_prefix);

			//perform testing
			vector<double> fold_margin_list = Test(test_id_list, mData);
			//add to test_result_map
			assert(fold_margin_list.size()==test_id_list.size());
			for (unsigned i = 0; i < fold_margin_list.size(); i++) {
				unsigned id = test_id_list[i];
				//pack all the result fields sequentially in a vector
				double target = test_target_list[i];
				double margin = fold_margin_list[i];
				double prediction = margin > 0 ? 1 : -1;
				vector<double> res;
				res.push_back(target);
				res.push_back(prediction);
				res.push_back(margin);
				//memoize the result vector with the test id
				oTestResultMap[id] = res;
			}
			cout << "Fold phase concluded in:" << endl;
		}

		vector<double> cv_prediction_list;
		vector<double> cv_target_list;
		string ofs_name = mParameters.mInputDataFileName + ".cv_predictions" + mParameters.mSuffix;
		ofstream ofs(ofs_name.c_str());
		//for all test ids read in order
		for (map<unsigned, vector<double> >::iterator it = oTestResultMap.begin(); it != oTestResultMap.end(); ++it) {
			unsigned id = it->first;
			//unpack the result fields from the result vector memoized with the test id
			double target = it->second[0];
			double prediction = it->second[1];
			double margin = it->second[2];
			ofs << id << " " << target << " " << prediction << " " << margin << endl;
			cv_prediction_list.push_back(prediction);
			cv_target_list.push_back(target);
		}
		cout << SEP << endl << "Performance on data set in cross validation:" << endl;
		double bfmeasure = mCoreModel.OutputPerformanceMeasures(cout, cv_prediction_list, cv_target_list);

		cout << endl << "Instance id, true target, prediction and margin saved in file: " << ofs_name << endl;
		cout << "Crossvalidation concluded in:" << endl;
		return bfmeasure;
	}

	void LearningCurve() {
		srand(mParameters.mRandomSeed);
		cout << SEP << endl;
		cout << "Computing learning curve with " << mParameters.mLearningCurveNumPoints << " folds." << endl;
		ProgressBar pbt;
		pbt.Count();
		vector<pair<double, double> > prediction_list;
		vector<pair<double, double> > margin_list;

		unsigned size = mData.Size();
		//randomly shuffle indices
		vector<unsigned> data_id_list;
		for (unsigned i = 0; i < size; ++i)
			data_id_list.push_back(i);
		for (unsigned i = 0; i < size; ++i) {
			unsigned j = (double) rand() / RAND_MAX * size;
			swap(data_id_list[i], data_id_list[j]);
		}

		//build train-test split in the learning curve way
		vector<unsigned> test_id_list;
		//test data is the first fold
		for (unsigned i = 0; i < size / mParameters.mLearningCurveNumPoints; ++i) {
			unsigned id = data_id_list[i];
			test_id_list.push_back(id);
		}
		//sort the indices in order to guarantee sequential file access
		sort(test_id_list.begin(), test_id_list.end());
		//extract target list for testing
		vector<double> test_target_list;
		for (unsigned i = 0; i < test_id_list.size(); i++) {
			unsigned id = test_id_list[i];
			test_target_list.push_back(mData.mTargetList[id]);
		}
		//training data is built incrementally adding 1/mParameters.mLearningCurveNumPoints * size instances
		for (unsigned f = 1; f < mParameters.mLearningCurveNumPoints; f++) { //NOTE: start from 1 as the first fold is used for the test data
			ProgressBar pbf;
			pbf.Count();
			cout << SEP << endl;
			cout << TAB << TAB << "Fold: " << f << " of " << mParameters.mLearningCurveNumPoints - 1 << endl;
			//generate the training set
			vector<unsigned> train_id_list;
			for (unsigned i = size / mParameters.mLearningCurveNumPoints; i < size * (f + 1) / mParameters.mLearningCurveNumPoints; ++i) {
				unsigned id = data_id_list[i];
				train_id_list.push_back(id);
			}
			//sort the indices in order to guarantee sequential file access
			sort(train_id_list.begin(), train_id_list.end());
			//extract target list for training
			vector<double> train_target_list;
			for (unsigned i = 0; i < train_id_list.size(); i++) {
				unsigned id = train_id_list[i];
				train_target_list.push_back(mData.mTargetList[id]);
			}

			//perform training
			Train(train_target_list, train_id_list, mData);

			//compute predictions on test
			{
				vector<double> fold_margin_list = Test(test_id_list, mData);
				assert(fold_margin_list.size()==test_id_list.size());
				cout << SEP << endl << "Performance on test set:" << endl;
				mCoreModel.OutputPerformanceMeasures(cout, fold_margin_list, test_target_list);

				//save results to file
				string ofs_name = mParameters.mInputDataFileName + ".lc_predictions_test_fold_" + stream_cast<string>(f) + mParameters.mSuffix;
				ofstream ofs(ofs_name.c_str());
				for (unsigned i = 0; i < fold_margin_list.size(); i++) {
					unsigned id = test_id_list[i];
					double target = test_target_list[i];
					double margin = fold_margin_list[i];
					double prediction = margin > 0 ? 1 : -1;
					ofs << id << " " << target << " " << prediction << " " << margin << endl;
				}
				ofs.close();
				cout << endl << "Instance id, true target, prediction and margin saved in file: " << ofs_name << endl;
			}
			//compute predictions on train
			{
				vector<double> fold_margin_list = Test(train_id_list, mData);
				assert(fold_margin_list.size()==train_id_list.size());
				cout << SEP << endl << "Performance on train set:" << endl;
				mCoreModel.OutputPerformanceMeasures(cout, fold_margin_list, train_target_list);

				//save results to file
				string ofs_name = mParameters.mInputDataFileName + ".lc_predictions_train_fold_" + stream_cast<string>(f) + mParameters.mSuffix;
				ofstream ofs(ofs_name.c_str());
				for (unsigned i = 0; i < fold_margin_list.size(); i++) {
					unsigned id = train_id_list[i];
					double target = train_target_list[i];
					double margin = fold_margin_list[i];
					double prediction = margin > 0 ? 1 : -1;
					ofs << id << " " << target << " " << prediction << " " << margin << endl;
				}
				ofs.close();
				cout << endl << "Instance id, true target, prediction and margin saved in file: " << ofs_name << endl;
			}

			cout << "Fold phase concluded in:" << endl;
		}
		cout << "Learning curve concluded in:" << endl;
	}

	void Matrix() {
		cout << SEP << endl << "Gram matrix phase" << endl << SEP << endl;
		ProgressBar pb;
		pb.Count();
		string output_filename = mParameters.mInputDataFileName + ".mtx" + mParameters.mSuffix;
		ofstream ofs;
		ofs.open(output_filename.c_str());
		if (!ofs) throw range_error("ERROR2.16: Cannot open file:" + output_filename);
		{
			cout << "Computing Gram matrix for " << mData.Size() << " instances." << endl;
			ProgressBar ppb;
			for (unsigned i = 0; i < mData.Size(); ++i) {
				for (unsigned j = 0; j < mData.Size(); ++j)
					ofs << mData.ComputeKernel(i, j) << " ";
				ofs << endl;
				ppb.Count();
			}
		}
		cout << "Gram matrix saved in file " << output_filename << endl;
	}

	void ComputeLocalMultiDimensionalScaling() {
		cout << SEP << endl << "Local Multi Dimensional Scaling phase" << endl << SEP << endl;

		LocalMultiDimensionalScaling M;
		M.Init(&mParameters, &mData);
		vector<FVector> x_list;
		M.ComputeLocalMultiDimensionalScaling(x_list);
		assert(x_list.size()==mData.Size());
		M.SaveEmbedding(x_list);
		M.SaveDistortion(x_list);
		M.SaveNeighbourhoodList();
	}

	void SaveModel(string aLocalSuffix = string()) {
		string filename = mParameters.mModelFileName + aLocalSuffix + mParameters.mSuffix;
		ofstream ofs;
		ofs.open(filename.c_str());
		if (!ofs) throw range_error("ERROR2.22: Cannot open file:" + filename);
		cout << endl << "Saving model file: " << filename << endl;
		mCoreModel.Save(ofs);
	}

	void LoadModel() {
		string filename = mParameters.mModelFileName + mParameters.mSuffix;
		ifstream ifs;
		ifs.open(filename.c_str());
		if (!ifs) throw range_error("ERROR2.23: Cannot open file:" + filename);
		cout << endl << "Loading model file: " << filename << endl;
		mCoreModel.Load(ifs);
		mCoreModel.OutputModelInfo();
		cout << endl;
	}

	void Train(vector<double> aTargetList, vector<unsigned> aTrainSetIDList, Data& aData) {
		assert(aTargetList.size()==aTrainSetIDList.size());
		//wrapper for semi-supervised case: self-training
		//assume unsupervised material receives 0 target
		//filter the unsupervised material and put it in separate lists
		//iterate:
		//train on supervised and test on unsupervised
		//replace 0 target with prediction
		vector<double> target_list(aTargetList);
		vector<unsigned> train_supervised_id_list;
		vector<double> train_supervised_target_list;
		vector<unsigned> train_unsupervised_id_list;
		for (unsigned i = 0; i < target_list.size(); ++i) {
			unsigned id = aTrainSetIDList[i];
			double target = target_list[i];
			if (target != 0) {
				train_supervised_id_list.push_back(id);
				train_supervised_target_list.push_back(target);
			} else {
				train_unsupervised_id_list.push_back(id);
			}
		}
		if (train_unsupervised_id_list.size() > 0) {
			//if unsupervised material is present then
			//train on supervised material
			cout << endl << "Semisupervised training on " << target_list.size() << " instances" << endl;
			cout << TAB << "supervised instances: " << train_supervised_id_list.size() << " (" << 100 * train_supervised_id_list.size() / (double) target_list.size() << "%)" << endl;
			cout << TAB << "unsupervised instances: " << train_unsupervised_id_list.size() << " (" << 100 * train_unsupervised_id_list.size() / (double) target_list.size() << "%)" << endl;
			vector<double> margin_list = mCoreModel.Train(train_supervised_target_list, train_supervised_id_list, aData);

			//repeat for a predefined number of iteration
			for (unsigned iteration = 0; iteration < mParameters.mSemiSupervisedNumIterations; iteration++) {
				//test on unsupervised material
				cout << endl << TAB << "Iteration " << iteration + 1 << "/" << mParameters.mSemiSupervisedNumIterations << endl;
				cout << "Testing on unsupervised instances: " << train_unsupervised_id_list.size() << endl;
				vector<double> margin_list = Test(train_unsupervised_id_list, aData);
				//find high and low threshold for margin (i.e. high confidence predictions)
				vector<double> sorted_margin_list;
				vector<double> sorted_positive_margin_list;
				vector<double> sorted_negative_margin_list;
				for (unsigned i = 0; i < margin_list.size(); ++i) {
					sorted_margin_list.push_back(margin_list[i]);
					if (margin_list[i] > 0) sorted_positive_margin_list.push_back(margin_list[i]);
					else sorted_negative_margin_list.push_back(margin_list[i]);
				}

				unsigned high_threshold_id, low_threshold_id;
				double high_threshold, low_threshold;
				cout
						<< "Predicted class distribution:  +:"
						<< sorted_positive_margin_list.size()
						<< " ("
						<< 100 * sorted_positive_margin_list.size() / (double) train_unsupervised_id_list.size()
						<< " %)"
						<< " -:"
						<< sorted_negative_margin_list.size()
						<< " ("
						<< 100 * sorted_negative_margin_list.size() / (double) train_unsupervised_id_list.size()
						<< " %)"
						<< endl;
				if (sorted_positive_margin_list.size() == 0 || sorted_negative_margin_list.size() == 0) {
					cout << "Warning: margins are one sided. Proceeding to use margin rank irrespectively of margin sign. Retaining " << mParameters.mSemiSupervisedThreshold / 2 * 100 << "% of most reliable predictions" << endl;
					sort(sorted_margin_list.begin(), sorted_margin_list.end());
					high_threshold_id = (sorted_margin_list.size() - 1) * (1 - mParameters.mSemiSupervisedThreshold / 2);
					low_threshold_id = (sorted_margin_list.size() - 1) * (mParameters.mSemiSupervisedThreshold / 2);
					high_threshold = sorted_margin_list.size() > 0 ? sorted_margin_list[high_threshold_id] : sorted_margin_list[sorted_margin_list.size() - 1];
					low_threshold = sorted_margin_list.size() > 0 ? sorted_margin_list[low_threshold_id] : sorted_margin_list[0];
				} else {
					sort(sorted_positive_margin_list.begin(), sorted_positive_margin_list.end());
					sort(sorted_negative_margin_list.begin(), sorted_negative_margin_list.end());
					high_threshold_id = (sorted_positive_margin_list.size() - 1) * (1 - mParameters.mSemiSupervisedThreshold);
					low_threshold_id = (sorted_negative_margin_list.size() - 1) * (mParameters.mSemiSupervisedThreshold);
					high_threshold = sorted_positive_margin_list.size() > 0 ? sorted_positive_margin_list[high_threshold_id] : sorted_negative_margin_list[sorted_negative_margin_list.size() - 1];
					low_threshold = sorted_negative_margin_list.size() > 0 ? sorted_negative_margin_list[low_threshold_id] : sorted_positive_margin_list[0];
				}
				cout << "Low score threshold:" << low_threshold << " High score threshold:" << high_threshold << endl;
				//replace 0 target with predicted target only for high confidence predictions
				map<unsigned, double> semi_supervise_augmented_target_map;
				for (unsigned i = 0; i < target_list.size(); ++i) //copy target for supervised instances
					if (target_list[i] != 0) {
						unsigned id = aTrainSetIDList[i];
						semi_supervise_augmented_target_map[id] = target_list[i];
					}
				unsigned counter_p = 0;
				unsigned counter_n = 0;
				for (unsigned i = 0; i < train_unsupervised_id_list.size(); ++i) { //copy prediction for unsupervised instances
					unsigned id = train_unsupervised_id_list[i];
					double margin = margin_list[i];
					double predicted_target;
					bool margin_test = false;
					if (mParameters.mSemiSupervisedInduceOnlyPositive) {
						margin_test = (margin >= high_threshold);
						predicted_target = 1;
					} else if (mParameters.mSemiSupervisedInduceOnlyNegative) {
						margin_test = (margin <= low_threshold);
						predicted_target = -1;
					} else {
						margin_test = (margin <= low_threshold || margin >= high_threshold);
						predicted_target = margin <= low_threshold ? -1 : 1;
					}
					if (margin_test == true) {
						assert( semi_supervise_augmented_target_map.count(id)==0);
						semi_supervise_augmented_target_map[id] = predicted_target;
						if (predicted_target > 0) counter_p++;
						else counter_n++;
					}
				}
				if (mParameters.mSemiSupervisedInduceOnlyPositive) cout << "Adding only predicted positives" << endl;
				if (mParameters.mSemiSupervisedInduceOnlyNegative) cout << "Adding only predicted negatives" << endl;
				cout << "Added +:" << counter_p << " and -:" << counter_n << " instances from unsupervised set to training set of size " << train_supervised_id_list.size() << endl;
				//compose indices vectors for training instances and target
				vector<unsigned> train_semi_supervise_augmented_id_list;
				vector<double> train_semi_supervise_augmented_target_list;
				for (map<unsigned, double>::iterator it = semi_supervise_augmented_target_map.begin(); it != semi_supervise_augmented_target_map.end(); ++it) {
					unsigned id = it->first;
					double target = it->second;
					train_semi_supervise_augmented_id_list.push_back(id);
					train_semi_supervise_augmented_target_list.push_back(target);
				}
				//retrain
				margin_list = mCoreModel.Train(train_semi_supervise_augmented_target_list, train_semi_supervise_augmented_id_list, aData);
			}

		} else { //if no unsupervised material is present then train directly
			vector<double> margin_list = mCoreModel.Train(target_list, aTrainSetIDList, aData);
		}
	}

	vector<double> Test(vector<unsigned> aTestSetIDList, Data& aData) {
		vector<unsigned> testset_id_list = aTestSetIDList;
		vector<double> margin_list = mCoreModel.Test(testset_id_list, aData);
		return margin_list;
	}

};

//------------------------------------------------------------------------------------------------------------------------
int main(int argc, const char **argv) {
	cout << SEP << endl << PROG_CREDIT << endl << SEP << endl;
	try {
		Model model;
		model.Init(argc, argv);
		model.Exec();
	} catch (exception& e) {
		cerr << e.what() << endl;
	}
	return 0;
}