#ifndef __MEX_HELPER_H__

#ifndef NO_CTL
#include <ctl.h>
#else
#include <iostream>
#include <vector>
#endif
#include <mex.h>

void initEnv();
bool mexTypeCheck (const mxArray *mval, const char *wantType, int i=1, int j=1);

#if __GNUC__ > 3
template<class T> class myVector : public std::vector<double>
#else
template<class T> class myVector : public std::vector<T>
#endif
{
	public:
	ctl::oStream &write(ctl::oStream &os) const 
	{
		os<<(const std::vector<T> &)*this;
	} 
	
	ctl::iStream &read(ctl::iStream &is)
	{
		is >> (std::vector<T> &)*this;
	}

	myVector () : std::vector<T>()
	{
	}

	myVector (int n) : std::vector<T>(n) 
	{ 
	}

	myVector (std::vector<T> v) : std::vector<T>(v)
	{
	}

	const myVector<T> get (int i, int j) const
	{
		myVector<T> ret(j-i);
		for (int k=i;k<j;k++)
			ret[k-i] = (*this)[k];
		return ret;
	}

	const myVector<T> operator* (const T j) const
	{
		myVector<T> ret(size());
		for (int i=0;i<size();i++)
			ret[i] = (*this)[i]*j;
		return ret;
	}

	// Warning: v is always treated as a column vector
	const T operator* (const myVector<T> &v) const
	{
		if (size() != v.size())
			return -1;

		T ret = 0;
		for (int i=0;i<size();i++)
			ret += (*this)[i] * v[i];
		return ret;
	}

	const myVector<T> operator- (const myVector<T> &v) const
	{
		if (size() != v.size())
			return myVector<T>();

		myVector<T> ret(size());
		for (int i=0;i<size();i++)
			ret[i] = (*this)[i] - v[i];
		return ret;
	}
};

template<class T> std::ostream& operator<<(std::ostream& os, const myVector<T> &v)
{
	for (int i=0;i<v.size();i++)
	{
			os << "\t";
			if (v[i] >= 0)
				os << " ";
			os << v[i];
			if (i!=v.size()-1)
				os << "\n";
	}
	return os;
}

template<class T> class matrix
{
#ifdef NO_CTL
	unsigned long long n;
#else
	ctl::uint8 n;
#endif
	std::vector<T> data;

	public: 
#ifndef NO_CTL
	CTL_Template(matrix, (T), 1, tupel, (n, data), 2);
#endif

	const int N () const
	{
		return n;
	}

	const int M () const
	{
		return data.size()/n;
	}

	const int size () const
	{
		return data.size();
	}

	matrix () : n(0)
	{
	}

	matrix (int n, int m)
	{
		resize(n, m);
	}

	matrix (int n, int m, const T* values)
	{
		resize(n, m);
		for (int k=0;k<n*m;k++)
			data[k] = values[k];
	}

	void resize (int n, int m)
	{
		this->n = n;
		data.resize(n*m);
	}

	const T &get (int i, int j) const
	{
		return data[i+n*j];
	}

	const myVector<T> get (int i) const
	{
		myVector<T> ret(n);
		for (int j=0;j<n;j++)
			ret[j] = data[i+n*j];
		return ret;
	}

	const myVector<T> get (int i, int j, int j2) const
	{
		myVector<T> ret(j2-j);
		for (int k=j;k<j2;k++)
			ret[k-j] = get(i, k);
		return ret;
	}

	bool set (int i, int j, T &dat)
	{
		data[i+n*j] = dat;
		return true;
	}
	
	bool set (int i, const myVector<T> &dat)
	{
		for (int j=0;j<n;j++)
			data[i+n*j] = dat[j];
		return true;
	}

	const T &operator[] (int i) const
	{
		return data[i];
	}

	T &operator[] (int i)
	{
		return data[i];
	}

	const matrix<T> operator* (const matrix<T> &B) const
	{
		matrix<T> ret;
		
		if (this->N() != B.M())
		{
			std::cerr << "Matrix dimensions do not agree.\n";
			return ret;
		}

		ret.resize(B.N(), this->M());
		for (int i=0;i<B.N();i++)
			for (int j=0;j<this->M();j++)
			{
				T tmp = 0;
				for (int k=0;k<this->N();k++)
					tmp += this->get(k, j)*B.get(i, k);
				ret.set(i, j, tmp);
			}

		return ret;
	}

	void print (std::ostream &os) const
	{
		for (int i=0;i<data.size();i++)
		{
			if (data[i] >= 0)
				os << " ";
			os << data[i] << " ";
			if (i<data.size()-1 && (i+1) % n == 0)
				os << "\n";
		}
	}
};

template<class T> std::ostream& operator<<(std::ostream& os, const matrix<T> &m)
{
	m.print(os);
	return os;
}

// Compile-time assertion for unsupported conversions
template<bool> struct Assert;
template<> struct Assert<true> {};

template<class T> inline bool mexConvert (const T &cval, mxArray *&mval)
{
	Assert<false> Conversion_not_supported;
	return false;
}

template<> inline bool mexConvert (const std::vector<double> &vec, mxArray *&mval)
{
	mval = mxCreateDoubleMatrix(vec.size(), 1, mxREAL);
	double *values = mxGetPr(mval);
	for (int i=0;i<vec.size();i++)
		values[i] = vec[i];
	return true;
}

template<> inline bool mexConvert (const matrix<double> &mat, mxArray *&mval)
{
	mval = mxCreateDoubleMatrix(mat.N(), mat.M(), mxREAL);
	double *values = mxGetPr(mval);
	for (int i=0;i<mat.size();i++)
			values[i] = mat[i];
	return true;
}

template<> inline bool mexConvert (const std::vector<std::vector<double> > &vec,
	mxArray *&mval)
{
	if (vec.size() == 0)
		mval = mxCreateDoubleMatrix(0, 0, mxREAL);
	else
	{
		if ((vec.size() >= 2) && (vec[0].size() != vec[1].size()))
		{
			std::cerr << "Conversion error: input is not a matrix.\n";
			return false;
		}
		mval = mxCreateDoubleMatrix(vec[0].size(), vec.size(), mxREAL);
		double *values = mxGetPr(mval);
		for (int j=0;j<vec[0].size();j++)
			for (int i=0;i<vec.size();i++)
				values[i*vec[0].size()+j] = vec[i][j];
	}
	return true;
}

template<> inline bool mexConvert (const float &f, mxArray *&mval)
{
	mval = mxCreateDoubleScalar(f);
	return true;
}

template<> inline bool mexConvert (const bool &b, mxArray *&mval)
{
	mval = mxCreateDoubleScalar(b);
	return true;
}

template<> inline bool mexConvert (const double &d, mxArray *&mval)
{
	mval = mxCreateDoubleScalar(d);
	return true;
}

template<> inline bool mexConvert (const int &i, mxArray *&mval)
{
	mval = mxCreateDoubleScalar(i);
	return true;
}

template<> inline bool mexConvert (const std::string &str, mxArray *&mval)
{
	return false;
}

#ifndef NO_CTL
template<> inline bool mexConvert (const ctl::string &str, mxArray *&mval)
{
	return false;
}

template<class T> inline bool mexConvert (const ctl::vector<T> &vec, mxArray *&mval)
{
	std::cerr << "Error: dummy converter used at runtime.\n";
	return false;
}
#endif

template<class T> inline bool mexConvert (const mxArray *mval, T &cval)
{
	Assert<false> Conversion_not_supported;
	return false;
}

#ifndef NO_CTL
template<class T> inline bool mexConvert (const mxArray *mval, ctl::vector<T> &vec)
{
	std::cerr << "Error: dummy converter used at runtime.\n";
	return false;
}

template<class T> inline bool mexConvert (const mxArray *mval, ctl::reference<T> &ref)
{
	std::cerr << "Error: dummy converter used at runtime.\n";
	return false;
}
#endif

template<> inline bool mexConvert (const mxArray *mval, std::vector<double> &vec)
{
	int len = mxGetM(mval);
	vec.resize(len);
	double *values = mxGetPr(mval);
	for (int i=0;i<len;i++)
		vec[i] = values[i];
	return true;
}

template<> inline bool mexConvert (const mxArray *mval, myVector<double> &vec)
{
	std::vector<double> moo = std::vector<double>(vec);
	mexConvert(mval, moo);
	vec.resize(moo.size());
	for (int i=0;i<moo.size();i++)
		vec[i] = moo[i];
	return true;
}

template<> inline bool mexConvert (const myVector<double> &vec, mxArray *&mval)
{
	std::vector<double> moo = std::vector<double>(vec);
	return mexConvert(moo, mval);
}

template<> inline bool mexConvert (const mxArray *mval, matrix<double> &mat)
{
	int N = mxGetN(mval);
	int M = mxGetM(mval);
	mat.resize(N, M);
	double *values = mxGetPr(mval);
	for (int i=0;i<N*M;i++)
		mat[i] = values[i];
	return true;
}

template<> inline bool mexConvert (const mxArray *mval, std::vector<std::vector<double> > &vec)
{
	int N = mxGetN(mval);
	int M = mxGetM(mval);
	vec.resize(N);
	double *values = mxGetPr(mval);
	for (int i=0;i<N;i++)
	{
		vec.resize(M);
		for (int j=0;j<M;j++)
			vec[i][j] = values[i*M+j];
	}
	return true;
}

template<> bool inline mexConvert (const mxArray *mval, double &cval)
{
	if (!mexTypeCheck(mval, "double"))
		return false;
	cval = mxGetScalar(mval);
	return true;
}

template<> bool inline mexConvert (const mxArray *mval, float &cval)
{
	if (!mexTypeCheck(mval, "double"))
		return false;
	cval = mxGetScalar(mval);
	return true;
}

template<> bool inline mexConvert (const mxArray *mval, int &cval)
{
	if (!mexTypeCheck(mval, "double"))
		return false;
	cval = (int)mxGetScalar(mval);
	return true;
}

template<> bool inline mexConvert (const mxArray *mval, std::string &cval)
{
	return false;
}

#endif