nr3matlab.h

matlab自回归马尔可夫转换模型仿真估计与预测

/* nr3matlab.h */
// version 0.8
// This file is a version of nr3.h with hooks that
// make it easy to write Matlab mex files, in particular
// ones that use NR3 routines.
// See http://www.nr.com/nr3_matlab.html

#include "mex.h"

#define _CHECKBOUNDS_ 1

// all the system #include's we'll ever need
#include <fstream>
#include <cmath>
#include <complex>
#include <iostream>
#include <iomanip>
#include <vector>
#include <limits>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fcntl.h>
#include <string.h>
#include <ctype.h>
#include <typeinfo.h>

using namespace std;

// macro-like inline functions

template<class T>
inline T SQR(const T a) {return a*a;}

template<class T>
inline const T &MAX(const T &a, const T &b)
        {return b > a ? (b) : (a);}

inline float MAX(const double &a, const float &b)
        {return b > a ? (b) : float(a);}

inline float MAX(const float &a, const double &b)
        {return b > a ? float(b) : (a);}

template<class T>
inline const T &MIN(const T &a, const T &b)
        {return b < a ? (b) : (a);}

inline float MIN(const double &a, const float &b)
        {return b < a ? (b) : float(a);}

inline float MIN(const float &a, const double &b)
        {return b < a ? float(b) : (a);}

template<class T>
inline T SIGN(const T &a, const T &b)
	{return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}

inline float SIGN(const float &a, const double &b)
	{return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}

inline float SIGN(const double &a, const float &b)
	{return (float)(b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));}

template<class T>
inline void SWAP(T &a, T &b)
	{T dum=a; a=b; b=dum;}

// exception handling when executing underneath Matlab

#ifdef _MSC_VER
#define throw(message) \
{char msg[1024]; sprintf_s(msg,1024,"%s in file %s at line %d\n", \
message,__FILE__,__LINE__); mexErrMsgTxt(msg);}
#else
#define throw(message) \
{char msg[1024]; sprintf(msg,"%s in file %s at line %d\n", \
message,__FILE__,__LINE__); mexErrMsgTxt(msg);}
#endif

// basic type names (put here so can use for Matlab stuff)

typedef int Int; // 32 bit integer
typedef unsigned int Uint;

#ifdef _MSC_VER
typedef __int64 Llong; // 64 bit integer
typedef unsigned __int64 Ullong;
#else
typedef long long int Llong; // 64 bit integer
typedef unsigned long long int Ullong;
#endif

typedef char Char; // 8 bit integer
typedef unsigned char Uchar;

typedef double Doub; // default floating type
typedef long double Ldoub;

typedef complex<double> Complex; // default complex type

typedef bool Bool;

// NaN: you should test by verifying that (NaN != NaN) is true (see nr3.h)
static const Doub NaN = numeric_limits<Doub>::quiet_NaN();

// get mxClassID of any type T

template <class T> inline mxClassID mxT() {return mxUNKNOWN_CLASS;}
template <> inline mxClassID mxT<Doub>() {return mxDOUBLE_CLASS;}
template <> inline mxClassID mxT<float>() {return mxSINGLE_CLASS;}
template <> inline mxClassID mxT<Int>() {return mxINT32_CLASS;}
template <> inline mxClassID mxT<Uint>() {return mxUINT32_CLASS;}
template <> inline mxClassID mxT<Char>() {return mxCHAR_CLASS;}
template <> inline mxClassID mxT<Uchar>() {return mxUINT8_CLASS;}
template <> inline mxClassID mxT<Llong>() {return mxINT64_CLASS;}
template <> inline mxClassID mxT<Ullong>() {return mxUINT64_CLASS;}
template <> inline mxClassID mxT<Bool>() {
	if (sizeof(Bool)==1) return mxLOGICAL_CLASS;
	else throw("bool and mxLOGICAL_CLASS have incompatible sizes");
}
inline mxClassID mxT(const mxArray *p) {return mxGetClassID(p);} 

// functions to map Matlab scalars

template <class T>
const T& mxScalar(const mxArray *prhs) {
	if (mxGetClassID(prhs) != mxT<T>())
		throw("attempt to assign scalar ref to wrong type");
	return *(T*)mxGetData(prhs);
}

template <class T>
T& mxScalar(mxArray* &plhs) {
	plhs = mxCreateNumericMatrix(1,1,mxT<T>(),mxREAL);
	return *(T*)mxGetData(plhs);
}

template <class T>
const T& mxScalar(const char *varname) {
	const mxArray* mxptr = mexGetVariablePtr("base",varname);
	if (mxptr == NULL) throw("attempt to get nonexistent variable");
	if (mxGetClassID(mxptr) != mxT<T>())
		throw("attempt to assign scalar ref to wrong type");
	return *(T*)mxGetData(mxptr);
}

template <class T>
void mxScalar(T val, const char *varname) {
	mxClassID mxclass = mxT<T>();
	if (mxclass == mxUNKNOWN_CLASS) throw("no corresponding Matlab type");
	mxArray* mxdum = mxCreateNumericMatrix(1,1,mxclass,mxREAL);
	*(T*)mxGetData(mxdum) = val;
	if (mexPutVariable("base",varname,mxdum))
		throw("failed to send data to Matlab variable by name");
	mxDestroyArray(mxdum);
}



// Vector and Matrix Classes

template <class T>
class NRvector {
private:
	int nn;	// size of array. upper index is nn-1
	int own; // 1 if own data, 0 if Matlab's
	T *v;
public:
	NRvector();
	explicit NRvector(int n);		// Zero-based array
	NRvector(int n, const T &a);	//initialize to constant value
	NRvector(int n, const T *a);	// Initialize to array
	NRvector(const NRvector &rhs);	// Copy constructor
	NRvector & operator=(const NRvector &rhs);	//assignment
	typedef T value_type; // make T available externally
	inline T & operator[](const int i);	//i'th element
	inline const T & operator[](const int i) const;
	inline int size() const;
	void resize(int newn); // resize (contents not preserved)
	void assign(int newn, const T &a); // resize and assign a constant value
	~NRvector();

	NRvector(const mxArray *prhs); // map Matlab rhs to vector (read-only)
	NRvector(int n, mxArray* &plhs); // create Matlab lhs and map to vector
	NRvector(const char *varname); // import Matlab variable by name (read-only)	
	void put(const char *varname); // export vector to a named Matlab variable

};

template <class T>
NRvector<T>::NRvector(const mxArray* prhs) : own(0) {
	if (mxGetClassID(prhs) != mxT<T>())
		throw("constructing VecDoub from a different Matlab type prhs");
	nn = mxGetNumberOfElements(prhs);
	v = (T*)mxGetData(prhs);
}

template <class T>
NRvector<T>::NRvector(int n, mxArray* &plhs) : nn(n), own(0) {
	mxClassID mxclass = mxT<T>();
	if (mxclass == mxUNKNOWN_CLASS) throw("no corresponding Matlab type for plhs");
	plhs = mxCreateNumericMatrix(1,nn,mxclass,mxREAL);
	v = (T*)mxGetData(plhs);
}

template <class T>
NRvector<T>::NRvector(const char *varname) : own(0) {
	const mxArray* mxptr = mexGetVariablePtr("base",varname);
	if (mxptr == NULL) throw("attempt to get nonexistent variable");
	if (mxGetClassID(mxptr) != mxT<T>())
		throw("constructing a VecDoub from a different Matlab type");
	nn = mxGetNumberOfElements(mxptr);
	v = (T*)mxGetData(mxptr);
}

template <class T>
void NRvector<T>::put(const char *varname) {
	mxClassID mxclass = mxT<T>();
	if (mxclass == mxUNKNOWN_CLASS) throw("no corresponding Matlab type");
	mxArray* mxdum = mxCreateNumericMatrix(1,1,mxclass,mxREAL);
	void* sav = mxGetData(mxdum);
	mxSetN(mxdum,nn);
	mxSetData(mxdum,v);
	if (mexPutVariable("base",varname,mxdum))
		throw("failed to send data to Matlab variable by name");
	mxSetData(mxdum,sav);
	mxSetN(mxdum,1);
	mxDestroyArray(mxdum);
}

template <class T>
NRvector<T>::NRvector() : nn(0), own(1), v(NULL) {}

template <class T>
NRvector<T>::NRvector(int n) : nn(n), own(1), v(n>0 ? new T[n] : NULL) {}

template <class T>
NRvector<T>::NRvector(int n, const T& a) : nn(n), own(1), v(n>0 ? new T[n] : NULL)
{
	for(int i=0; i<n; i++) v[i] = a;
}

template <class T>
NRvector<T>::NRvector(int n, const T *a) : nn(n), own(1), v(n>0 ? new T[n] : NULL)
{
	for(int i=0; i<n; i++) v[i] = *a++;
}

template <class T>
NRvector<T>::NRvector(const NRvector<T> &rhs) : nn(rhs.nn), own(1), v(nn>0 ? new T[nn] : NULL)
{
	for(int i=0; i<nn; i++) v[i] = rhs[i];
}

template <class T>
NRvector<T> & NRvector<T>::operator=(const NRvector<T> &rhs)
// postcondition: normal assignment via copying has been performed;
//		if vector and rhs were different sizes, vector
//		has been resized to match the size of rhs
{
	if (this != &rhs)
	{
		if (nn != rhs.nn) {
			if (!own) throw("resize of mxArray by assignment not allowed");
			if (v != NULL) delete [] (v);
			nn=rhs.nn;
			v= nn>0 ? new T[nn] : NULL;
		}
		for (int i=0; i<nn; i++)
			v[i]=rhs[i];
	}
	return *this;
}

template <class T>
inline T & NRvector<T>::operator[](const int i)	//subscripting
{
#ifdef _CHECKBOUNDS_
if (i<0 || i>=nn) {
	throw("NRvector subscript out of bounds");
}
#endif
	return v[i];
}

template <class T>
inline const T & NRvector<T>::operator[](const int i) const	//subscripting
{
#ifdef _CHECKBOUNDS_
if (i<0 || i>=nn) {
	throw("NRvector subscript out of bounds");
}
#endif
	return v[i];
}

template <class T>
inline int NRvector<T>::size() const
{
	return nn;
}

template <class T>
void NRvector<T>::resize(int newn)
{
	if (newn != nn) {
		if (!own) throw("resize of mxArray not allowed");
		if (v != NULL) delete[] (v);
		nn = newn;
		v = nn > 0 ? new T[nn] : NULL;
	}
}

template <class T>
void NRvector<T>::assign(int newn, const T& a)
{
	if (newn != nn) {
		if (!own) throw("resize of mxArray by assign method not allowed");
		if (v != NULL) delete[] (v);
		nn = newn;
		v = nn > 0 ? new T[nn] : NULL;
	}
	for (int i=0;i<nn;i++) v[i] = a;
}

template <class T>
NRvector<T>::~NRvector()
{
	if (own && v != NULL) delete[] (v);
}

// end of NRvector definitions

template <class T>
class NRmatrix {
private:
	int nn;
	int mm;
	int own; // owned by self 1, vs Matlab 0
	T **v;
public:
	NRmatrix();
	NRmatrix(int n, int m);			// Zero-based array
	NRmatrix(int n, int m, const T &a);	//Initialize to constant
	NRmatrix(int n, int m, const T *a);	// Initialize to array
	NRmatrix(const NRmatrix &rhs);		// Copy constructor
	NRmatrix & operator=(const NRmatrix &rhs);	//assignment
	typedef T value_type; // make T available externally
	inline T* operator[](const int i);	//subscripting: pointer to row i
	inline const T* operator[](const int i) const;
	inline int nrows() const;