matrix.cpp

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

		for(i=0; i<nn; i++)
			x[i]=m2[i+kk][i];

	return x;
}


vector matrix::diag(int k)
{
	return ::diag(*this,k);
}


//returns the sum of every element in the same row of input matrix.

vector sum(matrix  & m)
{
	int nr=m.mrow();
    vector x(nr);
	for(int i=0; i<nr; i++)
		x[i]=sum(m[i]);
	return x;
}


vector matrix::sum()
{
	int nr=mrow();
    vector x(nr);
	for(int i=0; i<nr; i++)
		x[i]=::sum((*this)[i]);
	return x;
}


matrix matrix::trans()
{
	//allocate the new matrix<double>
	int len=mrow();
	int ven=mcol();
    matrix v(ven,len);
	//do the operations
	for(int i=0; i<ven; i++)
		for(int j=0; j<len; j++)
		   v[i][j]=(*this)[j][i];
	return v;
}

matrix trans(matrix & v1)
{
	//allocate the new matrix<double>
	int len=v1.p->length;
	int ven=vlen(v1.p->f[0]);
    matrix v(ven,len);
	//do the operations
	for(int i=0; i<ven; i++)
		for(int j=0; j<len; j++)
		   v[i][j]=v1[j][i];
	return v;
}


matrix matrix::abs()
{
	int mrows=p->length;
	int mcols=vlen(p->f[0]);
    matrix v(mrows,mcols);
	for(int i=0; i<mrows; i++)
		for(int j=0; j<mcols; j++)
			v[i][j]=::fabs((*this)[i][j]);
	return v;
}


matrix abs(matrix & x)
{
	int mrows=x.mrow();
	int mcols=x.mcol();
    matrix v(mrows,mcols);
	for(int i=0; i<mrows; i++)
		for(int j=0; j<mcols; j++)
			v[i][j]=fabs(x[i][j]);
	return v;
}

//returns a vector. if the i-th row has any non-zero element,the i-th element
//of return vector is 1. otherwise, it is 0.
vector any(matrix & v)
{
	int m=v.mrow();
	int n=v.mcol();
    vector x(m);
	for(int i=0; i<m; i++){
		x[i]=0;
		for(int j=0; j<n; j++){
			if(v[i][j] != 0){
				x[i]=1;
				break;
			}
		}
	}
	return x;
}


vector matrix::any()
{
	int m=mrow();
	int n=mcol();
    vector x(m);
	for(int i=0; i<m; i++){
		x[i]=0;
		for(int j=0; j<n; j++){
			if((*this)[i][j] != 0){
				x[i]=1;
				break;
			}
		}
	}
	return x;
}

vector mmax(matrix & v)
{
	int nr=v.mrow();
	int nc=v.mcol();
    vector xx(nr);
	for(int i=0; i<nr; i++){
		xx[i]=v[i][0];
		for(int j=1; j<nc; j++)
			if(v[i][j]>xx[i])
				xx[i]=v[i][j];
	}
	return xx;
}

vector matrix::mmax()
{
	int nr=mrow();
	int nc=mcol();
    vector xx(nr);
	for(int i=0; i<nr; i++){
		xx[i]=(*this)[i][0];
		for(int j=1; j<nc; j++)
			if((*this)[i][j]>xx[i])
				xx[i]=(*this)[i][j];
	}
	return xx;
}

//v=max(v1,v2): v[i][j]=max(v1[i][j],v2[i][j].
matrix mmax(matrix & v1,matrix & v2)
{
	int nr=v1.mrow();
	int nc=v1.mcol();
    if( (v2.mrow() !=nr) || (v2.mcol() !=nc) )
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    matrix xx(nr,nc);
	for(int i=0; i<nr; i++)
		for(int j=0; j<nc; j++){
			xx[i][j]=v1[i][j];
			if(v2[i][j]>v1[i][j])
				xx[i][j]=v2[i][j];
		}
	return xx;
}

matrix matrix::mmax(matrix & v2)
{
	int nr=mrow();
	int nc=mcol();
    if( (v2.mrow() !=nr) || (v2.mcol() !=nc) )
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    matrix xx(nr,nc);
	for(int i=0; i<nr; i++)
		for(int j=0; j<nc; j++){
			xx[i][j]=(*this)[i][j];
			if(v2[i][j]>(*this)[i][j])
				xx[i][j]=v2[i][j];
		}
	return xx;
}

vector mmin(matrix & v)
{
	int nr=v.mrow();
	int nc=v.mcol();
    vector xx(nr);
	for(int i=0; i<nr; i++){
		xx[i]=v[i][0];
		for(int j=1; j<nc; j++)
			if(v[i][j]<xx[i])
				xx[i]=v[i][j];
	}
	return xx;
}

vector matrix::mmin()
{
	int nr=mrow();
	int nc=mcol();
    vector xx(nr);
	for(int i=0; i<nr; i++){
		xx[i]=(*this)[i][0];
		for(int j=1; j<nc; j++)
			if((*this)[i][j]<xx[i])
				xx[i]=(*this)[i][j];
	}
	return xx;
}

//v=min(v1,v2): v[i][j]=min(v1[i][j],v2[i][j].
matrix mmin(matrix & v1,matrix & v2)
{
	int nr=v1.mrow();
	int nc=v1.mcol();
    if( (v2.mrow() !=nr) || (v2.mcol() !=nc) )
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    matrix xx(nr,nc);
	for(int i=0; i<nr; i++)
		for(int j=0; j<nc; j++){
			xx[i][j]=v1[i][j];
			if(v2[i][j]<v1[i][j])
				xx[i][j]=v2[i][j];
		}
	return xx;
}


matrix matrix::mmin(matrix & v2)
{
	int nr=mrow();
	int nc=mcol();
    if( (v2.mrow() !=nr) || (v2.mcol() !=nc) )
    {
        cout<<Merrorstring[EMAT_INVALIDORDER];
        throw Merrorstring[EMAT_INVALIDORDER];
    }
    matrix xx(nr,nc);
	for(int i=0; i<nr; i++)
		for(int j=0; j<nc; j++){
			xx[i][j]=(*this)[i][j];
			if(v2[i][j]<(*this)[i][j])
				xx[i][j]=v2[i][j];
		}
	return xx;
}


matrix eye(int n)
{
    matrix x(n,n);
	for(int i=0; i<n; i++)
		x[i][i]=1;
	return x;
}


matrix eye(int m,int n)
{
    matrix x(m,n);
	int mn=MIN(m,n);
	for(int i=0; i<mn; i++)
		x[i][i]=1;
	return x;
}


matrix eye(matrix & v)
{
	int m=v.mrow();
	int n=v.mcol();
    matrix x(m,n);
	int mn=MIN(m,n);
	for(int i=0; i<mn; i++)
		x[i][i]=1;
	return x;
}

matrix matrix::eye()
{
	int m=mrow();
	int n=mcol();
    matrix x(m,n);
	int mn=MIN(m,n);
	for(int i=0; i<mn; i++)
		x[i][i]=1;
	return x;
}


matrix &  matrix::rand()
{
	int nr=mrow();
	int nc=mcol();
	for(int i=0; i<nr; i++){
		for(int j=0; j<nc; j++)
			 (*this)[i][j]=(double)::rand()/(double)RAND_MAX;
	}
	return *this;
}


matrix &  rand(matrix & v)
{
	int nr=v.mrow();
	int nc=v.mcol();
	for(int i=0; i<nr; i++){
		for(int j=0; j<nc; j++)
			v[i][j]=(double)rand()/(double)RAND_MAX;
	}
	return v;
}


//constructs and returns a m*n random matrix.
matrix matrix::rand(int m,int n)
{
    matrix v(m,n);
	for(int i=0; i<m; i++){
		for(int j=0; j<n; j++)
			v[i][j]=(double)::rand()/(double)RAND_MAX;
	}
	return v;
}


/*
double det(matrix & mm)
{
	int numcol=mm.mrow();
    if(mm.mcol() !=numcol)
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
    }

    matrix work(numcol,numcol);
	int nn;
	double det;
	work=copy(mm);
	nn=Matdet(work,numcol,det);
	if(nn !=0){
        throw Merrorstring[(Matvec_errors)nn];
             ;
	}
	return det;
}


double matrix::det()
{
	int numcol=mrow();
    if(mcol() !=numcol)
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
    }

    matrix work(numcol,numcol);
	int nn;
	double det;
	work=copy();
	nn=Matdet(work,numcol,det);
	if(nn !=0){
        throw Merrorstring[(Matvec_errors)nn];
             ;
	}
	return det;
}


matrix matrix::inv()
{
	int numcol=mrow();
    if(mcol() !=numcol)
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
    }

    matrix work(numcol,numcol);
	double de;
	int i,j,nn,c1;
	long chkmax;
	chkmax=(long)numcol*numcol;
	if(chkmax*sizeof(double)>=64000l){
		work=copy();
		nn=Matinv(work,numcol);
	}
	else{
		size_t xxx=(size_t)numcol*numcol;
        double *inv=new double[xxx];
		for(i=0; i<numcol; i++){
			c1=i*numcol;
			for(j=0; j<numcol; j++)
				inv[c1+j]=(*this)[i][j];
		}
		nn=Matinv(inv,numcol);
		for(i=0; i<numcol; i++){
			c1=i*numcol;
			for(j=0; j<numcol; j++)
				work[i][j]=inv[c1+j];
		}
		delete inv;
	 }

	   //because Matinv() may provide a error result at some time
	   //(it reports sigular when the matrix is not sigular),
	   //therefore, we use a new method (consume more time) Matinv1()
	   //to assure the sigular result that comes from Matinv().
	  if((nn==EMAT_SINGULAR) || (nn==WEMAT_SINGULAR))
	  {
			 work=copy();
			 nn=Matinv1(work,numcol,de);
	   }
       if(nn !=0){
             if(nn>0)
        throw Merrorstring[(Matvec_errors)nn];
             else
        throw Mwarningstring[(Matvec_warnings)nn];
       }
   return work;
}


matrix inv(matrix & mm)
{
	 return mm.inv();
}


matrix inv(matrix & mm, double & det)
{
	   int numcol=mm.mrow();
       if(mm.mcol() !=numcol)
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
    }
       matrix work(numcol,numcol);
	   int nn;
	   work=copy(mm);
	   nn=Matinv1(work,numcol,det);
       if(nn !=0){
             if(nn>0)
        throw Merrorstring[(Matvec_errors)nn];
             else
        throw Mwarningstring[(Matvec_warnings)nn];
       }
	   return work;
}


matrix matrix::inv(double & det)
{
	   int numcol=mrow();
       if(mcol() !=numcol)
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
    }

       matrix work(numcol,numcol);
	   int nn;
	   work=copy();
	   nn=Matinv1(work,numcol,det);
       if(nn !=0){
             if(nn>0)
        throw Merrorstring[(Matvec_errors)nn];
             else
        throw Mwarningstring[(Matvec_warnings)nn];
       }
	   return work;
}

matrix inv1(matrix & mm)
{
	   int numcol=mm.mrow();
       if(mm.mcol() !=numcol)
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
    }

       matrix work(numcol,numcol);
	   double det;
	   int nn;
	   work=copy(mm);
	   nn=Matinv1(work,numcol,det);
       if(nn !=0){
             if(nn>0)
        throw Merrorstring[(Matvec_errors)nn];
             else
        throw Mwarningstring[(Matvec_warnings)nn];
       }
	   return work;
}
*/

void matrix::matmult(vector& u,vector& v)
{
    int m=u();
    int n=v();
    initial(m,n);
    double* f1=u.p->f;
    double* f2=v.p->f;
//    vector **f=p->f;

    for(int i=0;i<m;i++,f1++)
    {
        f2=v.p->f;
        for(int j=0;j<n;j++,f2++)
//        (*f)[j]=(*f1)*(*f2);
        (*p->f[i])[j]=(*f1)*(*f2);
    }
}
void matrix::rndnorm(double mu,double sigma)
{
	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 vector*[len];
		for(int i=0; i<len; i++)
            p->f[i]=new vector(ven);
		p->tmppointer=NULL;
	}
    vector **f=p->f;
	for(int i=0; i<len; i++)
        f[i]->rndnorm(mu,sigma);
}
/*
double norm(matrix& v,int n)
{
	int nr=v.mrow();
	int nc=v.mcol();
	double res;
    matrix yy(nc,nr);
    vector zz(nr),zp(nc);
	switch(n)
	{
		case 0:
			zz=sum(abs(v));
			res=max(zz);
			break;
		case 1:
			zp=sum(abs(trans(v)));
			res=max(zp);
			break;
		case 2:
			zz=svd(v);
			res=zz[0];
			break;
		default:
			if(n<0){
				yy=(trans(v))*(v);
				res=sqrt(sum(diag(yy)));
			}
			else
        throw Merrorstring[EMAT_NODEFINE];
	}

	return res;
}


double matrix::norm(int n)
{
	int nr=mrow();
	int nc=mcol();
	double res;
    matrix yy(nc,nr);
    vector zz(nr),zp(nc);
	switch(n)
	{
		case 0:
			zz=((abs()).sum());
			res=zz.max();
			break;
		case 1:
			zp=((trans()).abs()).sum();
			res=zp.max();
			break;
		case 2:
			zz=svd();
			res=zz[0];
			break;
		default:
			if(n<0){
				yy=(trans())*(*this);
				res=::sqrt((yy.diag()).sum());
			}
			else
            ;
        throw Merrorstring[EMAT_NODEFINE];
	}

	return res;
}
*/

vector minofrow(matrix & v)
{
    vector xx;
    xx=v.minofrow();
	return xx;
}

vector matrix::minofrow()
{
	int nr=mrow();
	int nc=mcol();
    vector xx(nr);
	for(int i=0; i<nr; i++){
		xx[i]=(*this)[i][0];
		for(int j=1; j<nc; j++)
			if((*this)[i][j]<xx[i])
				xx[i]=(*this)[i][j];
	}
	return xx;
}

vector minofcol(matrix & v)
{
    vector xx;
    xx=v.minofcol();
	return xx;
}

vector matrix::minofcol()
{
	int nr=mrow();
	int nc=mcol();
    vector xx(nc);
    for(int i=0; i<nc; i++){
        xx[i]=(*this)[0][i];
        for(int j=1; j<nr; j++)
            if((*this)[j][i]<xx[i])
                xx[i]=(*this)[j][i];
	}
	return xx;
}

vector maxofrow(matrix & v)
{
    vector xx;
    xx=v.maxofrow();
	return xx;
}

vector matrix::maxofrow()
{
	int nr=mrow();
	int nc=mcol();
    vector xx(nr);
	for(int i=0; i<nr; i++){
		xx[i]=(*this)[i][0];
		for(int j=1; j<nc; j++)
			if((*this)[i][j]>xx[i])
				xx[i]=(*this)[i][j];
	}
	return xx;
}

vector maxofcol(matrix & v)
{
    vector xx;
    xx=v.maxofcol();
	return xx;
}

vector matrix::maxofcol()
{
	int nr=mrow();
	int nc=mcol();
    vector xx(nc);
    for(int i=0; i<nc; i++){
        xx[i]=(*this)[0][i];
        for(int j=1; j<nr; j++)
            if((*this)[j][i]>xx[i])
                xx[i]=(*this)[j][i];
	}
	return xx;
}

int matrix::copy(matrix & m)
{
	int mr=m.mrow();
	int mc=m.mcol();
	int f=1;
//    int f=initial(mr,mc);
	for(int i=0; i<mr; i++)
	for(int j=0; j<mc; j++)
		 (*this)[i][j]=m(i,j);
	return f;
}