intmat.cpp

基于线性规划的回归支持向量机源程序

#include "Intmat.hpp"
//the realization of IntMatrix template class
//#include "materr.hpp"
/*
char *Merrorstring[ElastError]=
{
	"[[error in error reporting ????]]",
	"the dimension is negative",
    "the dimensions of two IntMatrix conflict",
	"out of data range",
    "the IntMatrix is singular",
	"iteration failure",
	"data type error",
	"the index out of define range",
	"the divisor is close to zero",
	"the method used has failed",
	"the data has no definition" ,
	"the condition is infinite",
	"the method has not been defined",
	"capacity exceeded in construction",
    "function need a square IntMatrix"
};

//WARNING !! the sequence of next strings is in reverse order of
//			 enum Marvec_warnings !!
char *Mwarningstring[-WlastWarning]=
{
	"[[Warning in warning reporting ????]]",
	"condition is too large,the result may be invaliable",
	"norm is clear to zero, result is inaccurate",
    "IntMatrix is near singular, inv() result is inaccurate"
};
*/
extern char **Mwarningstring;
extern char **Merrorstring;

IntMatrix::IntMatrix()
{
	  p=new mrep;
	  p->length=0;
	  p->f=NULL;
	  p->refcnt=1;
	  p->tmppointer=NULL;
}



IntMatrix::IntMatrix(int xsize, int ysize, int init)
{
      if((xsize<=0)||(ysize<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(int);
      if((xsize>=floor(UINT_MAX/4))||(ysize>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  p=new mrep;
	  p->length=xsize;
      p->f=new IntVector *[xsize];
	  for(int i=0; i<xsize; i++)
          p->f[i]=new IntVector(ysize,init);
	   p->refcnt=1;
	  p->tmppointer=NULL;
}

int IntMatrix::initial(int xsize, int ysize)
{
      if((xsize<=0)||(ysize<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(int);
      if((xsize>=floor(UINT_MAX/4))||(ysize>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  int retn=1;
    if(p->refcnt==1)
	{
        if( (mrow()==xsize)&&(mcol()==ysize) )
        return 1;
		if(p->f !=NULL)
		{
			int len=p->length;
			for(int i=0; i<len; i++)
				 delete p->f[i];
			if(p->tmppointer !=NULL)
				delete p->tmppointer;
			delete p->f;
		}
	}
    else
    {
        p->refcnt--;
        p=new mrep;
    }
    if(xsize && ysize)
    {
        p->length=xsize;
        p->f=new IntVector *[xsize];
        for(int i=0; i<xsize; i++)
        p->f[i]=new IntVector(ysize);
    }

    p->refcnt=1;
    p->tmppointer=NULL;
	return retn;
}

int IntMatrix::initrow(int xsize)
{
      if((xsize<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(int);
      if(xsize>=floor(UINT_MAX/4))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  int retn=1;
    if(p->refcnt==1)
	{
        if( mrow()==xsize)
        return 1;
		if(p->f !=NULL)
		{
			int len=p->length;
			for(int i=0; i<len; i++)
				 delete p->f[i];
			if(p->tmppointer !=NULL)
				delete p->tmppointer;
			delete p->f;
		}
	}
    else
    {
        p->refcnt--;
        p=new mrep;
    }
    p->length=xsize;
    if( xsize )
    {
        p->f=new IntVector *[xsize];
        for(int i=0; i<xsize; i++)
        p->f[i]=new IntVector;
    }

    p->refcnt=1;
    p->tmppointer=NULL;
	return retn;
}


IntMatrix::IntMatrix(int  **a,int xlength,int ylength)
{
      if((xlength<=0)||(ylength<=0))
      {
        cout<< Merrorstring[EMAT_INVALIDSIZE];
        throw Merrorstring[EMAT_INVALIDSIZE];
    }

      int tmp=UINT_MAX /sizeof(int);
      if((xlength>=floor(UINT_MAX/4))||(ylength>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  p=new mrep;
	  p->length=xlength;
      p->f=new IntVector *[xlength];
	  for(int i=0; i<xlength; i++)
          p->f[i]=new IntVector(a[i],ylength);
	   p->refcnt=1;
	  p->tmppointer=NULL;
}

//constructs a new IntMatrix with IntMatrix a1 and a2. the construction mode is
//controlled by variable select:
//select=0:           the IntMatrix A1(n1*m1) and B1(n2*m2) will be placed at
//                    main diagnal and the dimension of new IntMatrix is
//					  (n1+n2)*(m1+m2);
//select=1:           combines IntMatrix A1(n1*m1) and B1(n2*m2).
//                    the dimension of new IntMatrix is n*(m1+m2)
//					  (n=MAX(n1,n2));
//select=2(default):  combines IntMatrix A1(n1*m1) and B1(n2*m2).
//                    the dimension of new IntMatrix is (n1+n2)*m
//					  (m=MAX(m1,m2));
//other select:		  reports an error !

IntMatrix::IntMatrix(IntMatrix  & a1,
        IntMatrix  & a2,int select)
{
	  int xlength,ylength,minc,ninc,i,j;
	  int mrow1=a1.mrow();
	  int mrow2=a2.mrow();
	  int mcol1=a1.mcol();
	  int mcol2=a2.mcol();
	  switch(select)
	  {
		case 0:
			xlength=mrow1+mrow2;
			ylength=mcol1+mcol2;
			minc=mrow1;
			ninc=mcol1;
			break;
		case 1:
			xlength=MAX(mrow1,mrow2);
			ylength=mcol1+mcol2;
			minc=0;
			ninc=mcol1;
			break;
		case 2:
			xlength=mrow1+mrow2;
			ylength=MAX(mcol1,mcol2);
			minc=mrow1;
			ninc=0;
			break;
		default:
        cout<<Merrorstring[EMAT_ASSIGNDATAERR];
        throw Merrorstring[EMAT_ASSIGNDATAERR];
        ;
	  }
      int tmp=UINT_MAX /sizeof(int);
      if((xlength>=floor(UINT_MAX/4))||(ylength>=tmp))
      {
        cout<< Merrorstring[EMAT_EXCEEDCAPACITY]<<endl;
        throw Merrorstring[EMAT_EXCEEDCAPACITY];
    }
	  p=new mrep;
	  p->length=xlength;
      p->f=new IntVector *[p->length];
	  for(i=0; i<xlength; i++)
          p->f[i]=new IntVector(ylength);
	  for(i=0; i<mrow1; i++)
		  for(j=0; j<mcol1; j++)
			 p->f[i]->p->f[j]=a1[i][j];
	  for(i=0; i<mrow2; i++)
		  for(j=0; j<mcol2; j++)
			p->f[i+minc]->p->f[j+ninc]=a2[i][j];
	  p->refcnt=1;
	  p->tmppointer=NULL;
}

//constructs an subIntMatrix using IntMatrix a. the dimension of new IntMatrix is
//(xlast-xfirst)*(ylast-yfist). the first element is a[xfirst][yfirst].

IntMatrix:: IntMatrix(IntMatrix  & a,int xfirst,int xlast,
						int yfirst,int ylast)
{
	  int xsize=xlast-xfirst;
	  int ysize=ylast-yfirst;
      if((xsize<=0)||(ysize<=0))
      {
        cout<<Merrorstring[EMAT_ASSIGNDATAERR];
        throw Merrorstring[EMAT_ASSIGNDATAERR];
    }
      if((xfirst<0)||(yfirst<0)||(xlast>a.mrow())||(ylast>a.mcol()))
      {
        cout<< Merrorstring[EMAT_OUTOFRANGE];
        throw Merrorstring[EMAT_OUTOFRANGE];
    }

	  p=new mrep;
	  p->length=xsize;
      p->f=new IntVector *[xsize];
	  for(int i=0; i<xsize; i++)
          p->f[i]=new IntVector(a[xfirst+i],yfirst,ylast);
	   p->refcnt=1;
	  p->tmppointer=NULL;
}



IntMatrix::IntMatrix(IntMatrix  & x)
{
	x.p->refcnt++;
	p=x.p;
}



IntMatrix::~IntMatrix()
{
	if(--p->refcnt==0)
	{
		if(p->f !=NULL)
		{
			int len=p->length;
			for(int i=0; i<len; i++)
				 delete p->f[i];
			if(p->tmppointer !=NULL)
				delete p->tmppointer;
			delete p->f;
		}
		delete p;
	}
}


IntMatrix copy(IntMatrix  & m)
{
	int mr=m.mrow();
	int mc=m.mcol();
    IntMatrix x(mr,mc);
	for(int i=0; i<mr; i++)
         *(x.p->f[i])=copy(*(m.p->f[i]));
	return x;
}


IntMatrix IntMatrix::copy()
{
	int mr=mrow();
	int mc=mcol();
    IntMatrix x(mr,mc);
	for(int i=0; i<mr; i++)
		 *(x.p->f[i])=::copy(*(p->f[i]));
	return x;
}



IntMatrix::operator int  ** ()
{
	if(p->tmppointer==NULL)
	{
		int n=mrow();

        p->tmppointer=new int *[n];
		for(int i=0; i<n; i++)
			  p->tmppointer[i]=p->f[i]->p->f;
	}
	return p->tmppointer;
}



IntMatrix  &  IntMatrix::operator=(IntMatrix  & m)
{
	m.p->refcnt++;
	if(--p->refcnt==0)
	{
		int len=p->length;
		for(int i=0; i<len; i++)
			delete p->f[i];
		delete p->f;
		if(p->tmppointer !=NULL)
			delete p->tmppointer;
		delete p;
	}
	p=m.p;
	return *this;
}


IntMatrix  &  IntMatrix::operator=(IntVector  & x)
{
	int len=p->length;
	if(p->refcnt>1)
    {                    //x has been referenced by other IntMatrix.
		int ven=p->f[0]->p->length;
		p->refcnt--;
        p=new mrep;      //produce a new IntMatrix data structure.
		p->length=len;
		p->refcnt=1;
        p->f=new IntVector*[len];
		for(int i=0; i<len; i++)
            p->f[i]=new IntVector(ven);
		p->tmppointer=NULL;
	}
    IntVector **f=p->f;
	for(int i=0; i<len; i++)
		*f[i]=x;
	return *this;
}


IntMatrix  &  IntMatrix::operator=(int x)
{
	int len=p->length;
	if(p->refcnt>1)
	{
		int ven=p->f[0]->p->length;
		p->refcnt--;
		p=new mrep;
		p->length=len;
		p->refcnt=1;
        p->f=new IntVector*[len];
		for(int i=0; i<len; i++)
            p->f[i]=new IntVector(ven);
		p->tmppointer=NULL;
	}
    IntVector **f=p->f;
	for(int i=0; i<len; i++)
		*f[i]=x;
	return *this;
}


//m=m1+m2

IntMatrix operator+(IntMatrix  & m1,IntMatrix  & m2)
{
    //allocate the new IntMatrix
    if(m1!=m2)
    {
        cout<<"IntMatrix's dim are not equal";
    }
	int len=MIN(m1.p->length,m2.p->length);
	int ven=MIN(vlen(m1.p->f[0]),vlen(m2.p->f[0]));

    IntMatrix m(len,ven);
	//do the operations
    IntVector **f=m.p->f;
    IntVector **f1=m1.p->f;
    IntVector **f2=m2.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] + *f2[i];
	return m;
}

//m=m1-m2

IntMatrix operator-(IntMatrix  & m1,IntMatrix  & m2)
{
    //allocate the new IntMatrix
    if(m1!=m2)
    {
        cout<<"IntMatrix's dim are not equal";
    }
	int len=MIN(m1.p->length,m2.p->length);
	int ven=MIN(vlen(m1.p->f[0]),vlen(m2.p->f[0]));

    IntMatrix m(len,ven);
	//do the operations
    IntVector **f=m.p->f;
    IntVector **f1=m1.p->f;
    IntVector **f2=m2.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] - *f2[i];
	return m;
}


//m=m1*m2

IntMatrix operator*(IntMatrix  & m1,IntMatrix  & m2)
{
	//check the dimensions
	int i,j;
	int m1rows=m1.mrow();
	int m1cols=m1.mcol();
	int m2rows=m2.mrow();
	int m2cols=m2.mcol();
    if(m1cols != m2rows)
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    IntMatrix m(m1rows,m2cols);
    // create the IntMatrix with its row and columns flips
    IntMatrix flip(m2cols, m2rows);
	for(i=0; i<m2rows; i++)
		for(j=0; j<m2cols; j++)
			flip[j][i]=m2[i][j];

	//create the cross product
	for( i=0; i<m1rows; i++)
		for( j=0; j<m2cols; j++){
			m[i][j]=m1[i].dot(flip[j]);
		}
	return m;
}


//m=m1/m2

IntMatrix operator/(IntMatrix  & m1,IntMatrix  & m2)
{
    //allocate the new IntMatrix
    if(m1!=m2)
    {
        cout<<"IntMatrix's dim are not equal";
    }
	int len=MIN(m1.p->length,m2.p->length);
	int ven=MIN(vlen(m1.p->f[0]),vlen(m2.p->f[0]));

    IntMatrix m(len,ven);
	//do the operations
    IntVector **f=m.p->f;
    IntVector **f1=m1.p->f;
    IntVector **f2=m2.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] / *f2[i];
	return m;
}

//m=m1+x

IntMatrix operator+(IntMatrix  & m1,int x)
{
    //allocate the new IntMatrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    IntMatrix m(len,ven);
	//do the operations
    IntVector **f=m.p->f;
    IntVector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] + x;
	return m;
}

//m=m1-x

IntMatrix operator-(IntMatrix  & m1,int x)
{
    //allocate the new IntMatrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    IntMatrix m(len,ven);
	//do the operations
    IntVector **f=m.p->f;
    IntVector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] - x;
	return m;
}

//m=m1*x

IntMatrix operator*(IntMatrix  & m1,int x)
{
    //allocate the new IntMatrix
	int len=m1.p->length;
	int ven=vlen(m1.p->f[0]);

    IntMatrix m(len,ven);
	//do the operations
    IntVector **f=m.p->f;
    IntVector **f1=m1.p->f;
	for(int i=0; i<len; i++)
		*f[i]=*f1[i] *x;
	return m;
}


//v=m1*v2

IntVector operator*(IntMatrix  & m1,IntVector  & v2)
{
	//check the dimensions
	int m1rows=m1.mrow();
	int m1cols=m1.mcol();
	int m2rows=v2.vlen();
    if(m1cols != m2rows)
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    IntVector v(m1rows);

	//create the cross product
	for( int i=0; i<m1rows; i++){
			 v[i]=m1[i].dot(v2);
	}
	return v;
}