lmse.cpp

文件中包含多种模式识别常用的算法

#include<iostream.h>
#include<math.h>
#include<stdlib.h>

void mt(float *a,float *b,int row,int column)//矩阵求转置函数
{
	int i,j;//a是转置前矩阵 b是转置后的矩阵 row为矩阵的行数 column是矩阵的列数
    float temp;
    for(i=0;i<row;i++)
		for(j=0;j<column;j++)
		{
			temp=*(a+i*column+j);
			*(b+j*row+i)=temp;
		}
}

void subb(float *a,float *b,int n,float *c)//矩阵相减函数
{                                          //c中存放差矩阵
   int i;
   for(i=0;i<n;i++)
   {
	   *(c+i)=*(a+i)-*(b+i);
   }
}

void madd(float *a,float *b,int n,float *c)//矩阵相加函数
{                                          //c中存放和矩阵
   int i;
   for(i=0;i<n;i++)
   {
	   *(c+i)=*(a+i)+*(b+i);//fabs是浮点数绝对值函数
   }
}

void kmul(float *a,int n,int c,float *b)//矩阵数乘函数
{
   int i;
   for(i=0;i<n;i++)
   {
	   *(b+i)=(*(a+i))*c;
   }
}

void bmul(float *a,float *b,int m,int n,int k,float *result)//矩阵相乘矩阵
{ 
	int i,j,l,u;              //矩阵result=矩阵a*矩阵b
    for (i=0; i<=m-1; i++)    //m是矩阵a的行数，n是矩阵a的列数 k是矩阵b的列数
		for (j=0; j<=k-1; j++)
		{
			u=i*k+j;
			result[u]=0.0;
			for (l=0; l<=n-1; l++)
				result[u]=result[u]+a[i*n+l]*b[l*k+j];
		}
    return;
}

int inv(float *a,int n) //矩阵求逆函数 n是矩阵a的阶数
{
	int *is,*js,i,j,k,l,u,v;
    float d,p;
    
    is=new int [n];
    js=new int [n];
  
    for (k=0; k<=n-1; k++)
	{
		d=0.0;
        for (i=k; i<=n-1; i++)
			for (j=k; j<=n-1; j++)
			{
				l=i*n+j;
				p=fabs(a[l]);
				if (p>d) { d=p; is[k]=i; js[k]=j; }
			}
        if (d+1.0==1.0)
		{
			delete is;
			delete js;
			cout<<"error matrix can not inv"<<endl;
            return(0);
		}
        if (is[k]!=k)
			for (j=0; j<=n-1; j++)
			{
				u=k*n+j; v=is[k]*n+j;
				p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (js[k]!=k)
			for (i=0; i<=n-1; i++)
            {
				u=i*n+k; v=i*n+js[k];
				p=a[u]; a[u]=a[v]; a[v]=p;
            }
        l=k*n+k;
        a[l]=1.0/a[l];
        for (j=0; j<=n-1; j++)
		  if (j!=k)
		  { u=k*n+j; a[u]=a[u]*a[l]; }
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
              if (j!=k)
                {
				  u=i*n+j;
                  a[u]=a[u]-a[i*n+k]*a[k*n+j];
                }
        for (i=0; i<=n-1; i++)
          if (i!=k)
		  { u=i*n+k; a[u]=-a[u]*a[l];}
	}
    for (k=n-1; k>=0; k--)
	{
		if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            {
			  u=k*n+j; v=js[k]*n+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            {
			  u=i*n+k; v=i*n+is[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
	}
    delete is; delete js;
    return(1);
}


void main()
{
	float x[8][4]={{0,0,0,1},{1,0,0,1},{1,0,1,1},{1,1,0,1},{0,0,-1,-1},{0,-1,-1,-1},{0,-1,0,-1},{-1,-1,-1,-1}};
	float xt[4][8];
	float xm[4][4];
	float xn[4][8];
	float xw[8];
	float w[4];
	float wa[4];
	float b[8]={1,1,1,1,1,1,1,1};
	float ba[8];
	float e[8];
	float ef[8];
	int c=2;
	int i;
	int f0,f1,f2;
	
	//x的伪逆矩阵
	mt(&x[0][0],&xt[0][0],8,4);
	bmul(&xt[0][0],&x[0][0],4,8,4,&xm[0][0]);
	inv(&xm[0][0],4);
	bmul(&xm[0][0],&xt[0][0],4,4,8,&xn[0][0]);
	
	bmul(&xn[0][0],&b[0],4,8,1,&w[0]); //w的初值
	while(1)
	{
		bmul(&x[0][0],&w[0],8,4,1,&xw[0]);
		subb(&xw[0],&b[0],8,&e[0]); //e的值
		//根据e的值进行判别
		f0=0; f1=0; f2=0;
		for(i=0;i<8;i++)
		{
			if(e[i]==0) f0++;
			else if(e[i]>0) f1++;
			else f2++;
			ef[i]=fabs(e[i]);
		}
		if(f0==8) //e=0得解
			break;
		if(f1>0) //e>0继续迭代
		{
			bmul(&xn[0][0],&ef[0],4,8,1,&wa[0]);
			kmul(&wa[0],4,c,&wa[0]);
			madd(&w[0],&wa[0],4,&w[0]); //修正w

			madd(&e[0],&ef[0],8,&ba[0]);
			kmul(&ba[0],8,c,&ba[0]);
			madd(&b[0],&ba[0],8,&b[0]); //修正b
		}
		if(f2==8) //e的全部分量为负值，模式线性不可分
		{
			cout<<"pattern cannot be classified linearly!"<<endl;
			break;
		}
	}

	cout<<"the solution is:("<<w[0]<<","<<w[1]<<","<<w[2]<<","<<w[3]<<")"<<endl;
}