matlib.cpp

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

#include "Intmat.hpp"
#include "matrix.hpp"
#include "IntVect.hpp"

#include<math.h>
#include<conio.h>
#include "matlib.h"

extern char **Mwarningstring;
extern char **Merrorstring;

int SolveSysEqa(matrix& a,matrix& c)
{
    int i,j,k,k1,k2,k3;
    double p;
    if(! (a.mrow()&&a.mcol()&&c.mrow()&&c.mcol())  )
    return -1;
    if( a.mrow()!=a.mcol() )
    return -2;
    int n=a.mrow();
    int m=c.mcol();
    if(a.mrow()!=c.mrow())
    return -3;
    if (fabs(a(0,0))+1.0==1.0)
    {
        return 0;
    }
    for (i=1; i<=n-1; i++)
    {
        a(i,0)=a(i,0)/a(0,0);
    }
    for (i=1; i<=n-2; i++)
    {
//        u=i*n+i;
        for (j=1; j<=i; j++)
        {
//            v=i*n+j-1; l=(j-1)*n+j-1;
            a(i,i)=a(i,i)-a(i,j-1)*a(i,j-1)*a(j-1,j-1);
        }
        p=a(i,i);
        if (fabs(p)+1.0==1.0)
        {
            return 0;
        }
        for (k=i+1; k<=n-1; k++)
        {
//            u=k*n+i;
            for (j=1; j<=i; j++)
            {
//                v=k*n+j-1; l=i*n+j-1; w=(j-1)*n+j-1;
                a(k,i)=a(k,i)-a(k,j-1)*a(i,j-1)*a(j-1,j-1);
            }
            a(k,i)=a(k,i)/p;
        }
      }
//    u=n*n-1;
    for (j=1; j<=n-1; j++)
    {
//        v=(n-1)*n+j-1; w=(j-1)*n+j-1;
        a(n-1,n-1)=a(n-1,n-1)-a(n-1,j-1)*a(n-1,j-1)*a(j-1,j-1);
    }
    p=a(n-1,n-1);
    if(fabs(p)+1.0==1.0)
    return 0;
    for(j=0; j<=m-1; j++)
    for(i=1; i<=n-1; i++)
    {
//        u=i*m+j;
        for (k=1; k<=i; k++)
        {
//            v=i*n+k-1; w=(k-1)*m+j;
            c(i,j)=c(i,j)-a(i,k-1)*c(k-1,j);
        }
      }
    for(i=1; i<=n-1; i++)
    {
//        u=(i-1)*n+i-1;
        for(j=i; j<=n-1; j++)
        {
//            v=(i-1)*n+j; w=j*n+i-1;
            a(i-1,j)=a(i-1,i-1)*a(j,i-1);
        }
    }
    for(j=0; j<=m-1; j++)
    {
//        u=(n-1)*m+j;
        c(n-1,j)=c(n-1,j)/p;
        for (k=1; k<=n-1; k++)
        {
            k1=n-k; k3=k1-1;
//            u=k3*m+j;
            for (k2=k1; k2<=n-1; k2++)
            {
//                v=k3*n+k2; w=k2*m+j;
                c(k3,j)=c(k3,j)-a(k3,k2)*c(k2,j);
            }
            c(k3,j)=c(k3,j)/a(k3,k3);
          }
    }
    return 1;
}
double SolveEqa(matrix& a,matrix& c)
{
    int amr=a.mrow();
    int amc=a.mcol();
    int cmr=c.mrow();
    int cmc=c.mcol();
    if(! (amr&&amc&&cmr&&cmc)  )
    return -10;
    if(amr!=cmr)
    return -10;
    if(amc>amr)
    return -10;
    int mr=amr;
    int mc=amc+cmc;
    matrix b(mr,mc);
    for(int i=0;i<mr;i++)
    {
        for(int j=0;j<amc;j++)
        b(i,j)=a(i,j);
        for(j=0;j<cmc;j++)
        b(i,amc+j)=c(i,j);
    }
    qr1(b);
    c.initial(amc,cmc);
    for(i=amc-1;i>=0;i--)
    {
        for(int j=0;j<cmc;j++)
        {
            c(i,j)=b(i,j+amc)/b(i,i);
            for(int k=i+1;k<amc;k++)
            c(i,j)-=b(i,k)*c(k,j)/b(i,i);
        }
    }
    double err=0;
    mr=MIN(amr,cmc+amc);
    for(i=cmc;i<mr;i++)
    for(int j=i;j<mc;j++)
    err+=b(i,j)*b(i,j);
    return err;
}

//col(a)=vlen(v),solve eqation x*a=v,vlen(x)=mrow(a)
double SolveREqa(matrix& a,vector& v,vector& x)
{
    int m=mrow(a);
    int n=mcol(a);
    double delta=-1.0;
    if(n!=vlen(v))
    return delta;
    matrix b(n,m+1);
    for(int i=0;i<n;i++)
    for(int j=0;j<m;j++)
    b(i,j)=a(j,i);
    for(i=0;i<m;i++)
    b(i,m)=v[i];
    qr1(b);
    delta=b(m,m);
    x.initial(m);
    x[m-1]=b(m-1,m)/b(m-1,m-1);
    for(i=m-2;i>=0;i--)
    {
        x[i]=b(i,m);
        for(int j=i+1;j<m;j++)
        x[i]-=x[j]*b(i,j);
        x[i]/=b(i,i);
    }
    return delta;
}
//mrow(a)=vlen(v),a*x=v,vlen(x)=mcol(a)
double SolevLEqa(matrix& a,vector& v,vector& x)
{
    int m=mrow(a);
    int n=mcol(a);
    double delta=-1.0;
    if(m!=vlen(v))
    return delta;
    matrix b(m,n+1);
    for(int i=0;i<m;i++)
    for(int j=0;j<n;j++)
    b(i,j)=a(i,j);
    for(i=0;i<m;i++)
    b(i,n)=v[i];
    qr1(b);
    delta=b(n,n);
    x.initial(n);
    x[n-1]=b(n-1,n)/b(n-1,n-1);
    for(i=n-2;i>=0;i--)
    {
        x[i]=b(i,n);
        for(int j=i+1;j<n;j++)
        x[i]-=x[j]*b(i,j);
        x[i]/=b(i,i);
    }
    return delta;
}

int eastrq(matrix& a,matrix& q,vector& b,vector& c)
{
    int i,j,k;
    double h,f,g,h2;
    if(! (a.mrow()&&a.mcol()) )
    return -1;
    if(a.mrow()!=a.mcol() )
    return -2;
    int n=a.mrow();
    b.initial(n);
    c.initial(n);
    q=copy(a);
    for (i=n-1; i>=1; i--)
    {
        h=0.0;
        if (i>1)
        for (k=0; k<=i-1; k++)
        h=h+q(i,k)*q(i,k);
        if (h+1.0==1.0)
        {
            c[i]=0.0;
            if (i==1) c[i]=q(i,i-1);
            b[i]=0.0;
        }
        else
        {
            c[i]=sqrt(h);
            if (q(i,i-1)>0.0) c[i]=-c[i];
            h=h-q(i,i-1)*c[i];
            q(i,i-1)=q(i,i-1)-c[i];
            f=0.0;
            for (j=0; j<=i-1; j++)
            {
                q(j,i)=q(i,j)/h;
                g=0.0;
                for (k=0; k<=j; k++)
                g=g+q(j,k)*q(i,k);
                if (j+1<=i-1)
                for (k=j+1; k<=i-1; k++)
                g=g+q(k,j)*q(i,k);
                c[j]=g/h;
                f=f+g*q(j,i);
            }
            h2=f/(h+h);
            for (j=0; j<=i-1; j++)
            {
                f=q(i,j);
                g=c[j]-h2*f;
                c[j]=g;
                for (k=0; k<=j; k++)
                {
                    q(j,k)=q(j,k)-f*c[k]-g*q(i,k);
                }
            }
            b[i]=h;
        }
    }
    for(i=0;i<n-1;i++)
    c[i]=c[i+1];
    c[n-1]=0.0;
        b[0]=0.0;
        for (i=0; i<=n-1; i++)
        {
            if ((b[i]!=0.0)&&(i-1>=0))
            for (j=0; j<=i-1; j++)
            {
                g=0.0;
                for (k=0; k<=i-1; k++)
                g=g+q(i,k)*q(k,j);
                for (k=0; k<=i-1; k++)
                {
                    q(k,j)=q(k,j)-g*q(k,i);
                }
            }
            b[i]=q(i,i); q(i,i)=1.0;
            if (i-1>=0)
            for (j=0; j<=i-1; j++)
            { q(i,j)=0.0; q(j,i)=0.0;}
        }
        return 1;
}
int ebstq(vector b,vector c,matrix& q,double eps,int l)
{
    int i,j,k,m,it,n;
    double d,f,h,g,p,r,e,s;
    n=b();
    if(!n || (b()!=c()) )
    return 0;
    c[n-1]=0.0;
    d=0.0; f=0.0;
    for (j=0; j<=n-1; j++)
    {
        it=0;
        h=eps*(fabs(b[j])+fabs(c[j]));
        if (h>d) d=h;
        m=j;
        while ((m<=n-1)&&(fabs(c[m])>d)) m=m+1;
        if (m!=j)
          { do
              {
                if (it==l)
                return -1;
                it=it+1;
                g=b[j];
                p=(b[j+1]-g)/(2.0*c[j]);
                r=sqrt(p*p+1.0);
                if (p>=0.0) b[j]=c[j]/(p+r);
                else b[j]=c[j]/(p-r);
                h=g-b[j];
                for (i=j+1; i<=n-1; i++)
                  b[i]=b[i]-h;
                f=f+h; p=b[m]; e=1.0; s=0.0;
                for (i=m-1; i>=j; i--)
                  { g=e*c[i]; h=e*p;
                    if (fabs(p)>=fabs(c[i]))
                      { e=c[i]/p; r=sqrt(e*e+1.0);
                        c[i+1]=s*p*r; s=e/r; e=1.0/r;
                      }
                    else
		      { e=p/c[i]; r=sqrt(e*e+1.0);
                        c[i+1]=s*c[i]*r;
                        s=1.0/r; e=e/r;
                      }
                    p=e*b[i]-s*g;
                    b[i+1]=h+s*(e*g+s*b[i]);
                    for (k=0; k<=n-1; k++)
                    {
//                        u=k*n+i+1; v=u-1;
                        h=q(k,i+1); q(k,i+1)=s*q(k,i)+e*h;
                        q(k,i)=e*q(k,i)-s*h;
                      }
                  }
                c[j]=s*p; b[j]=e*p;
              }
            while (fabs(c[j])>d);
          }
        b[j]=b[j]+f;
      }
    for (i=0; i<=n-1; i++)
      { k=i; p=b[i];
        if (i+1<=n-1)
          { j=i+1;
            while ((j<=n-1)&&(b[j]<=p))
              { k=j; p=b[j]; j=j+1;}
          }
        if (k!=i)
          { b[k]=b[i]; b[i]=p;
            for (j=0; j<=n-1; j++)
              {
                //u=j*n+i; v=j*n+k;
                p=q(j,i); q(j,i)=q(j,k); q(j,k)=p;
              }
          }
      }
    return 1;
}

int SysEig(matrix& a, matrix& p, vector& d)
{
    vector c;
    double eps=0.0000001;
    int l=100000;
    if(a.mrow()!=a.mcol())
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
        return 0;
    }
    double dmin,d2;
    dmin=a(0,0);
    for(int i=0;i<a.mrow();i++)
    {
        d2=0;
        for(int j=0;j<a.mcol();j++)
        {
            if(i==j)
            d2+=fabs(a(i,j));
            else
            d2-=fabs(a(i,j));
        }
        if(dmin>d2)
        dmin=d2;
    }
    if(dmin<0.1)
    for(int i=0;i<a.mrow();i++)
    a(i,i)+=(0.1+fabs(dmin));
try
{
    if(eastrq(a,p,d,c)>0)
    {
	ebstq(d,c,p,eps,l);
    }
    if(dmin<0.1)
    for(int i=0;i<a.mrow();i++)
    d[i]-=(0.1+fabs(dmin));
}
catch(char* str)
{
    cout<<"SysEig="<<str<<endl;
}
catch(...)
{
    cout<<"SysEig other exception"<<endl;
}
    return 1;
}


int SelfCov(vector& x,vector& r,int n)
{
    double me=0.0;
    vector y;
    if(!x())
    return 0;
    y=copy(x);
    me=sum(x)/vlen(x);
    y-=me;
    r.initial(n);
    for(int i=0;i<n;i++)
    {
        r[i]=0.0;
        for(int j=0;j<=(x()-1-i);j++)
        r[i]=r[i]*j/(j+1)+y[j]*y[j+i]/(j+1);
        r[i]=r[i]*(x()-i)/x();
    }
    return 1;
}
int PartialCorrByLS(matrix& a,vector& pr)
{
    int n,m;
    m=a.mrow();
    n=a.mcol();
    if(m<n)
    return 0;
    pr.initial(n-1);
    for( int i=0;i<(n-1);i++)
    {
        if(fabs(a(i,i))<=0.00000001)
        return 0;
        pr[i]=a(i,n-1)/a(i,i);
    }
    return 1;
}

int qr(matrix& a)
{
    int i,j,k,nn,jj,m,n;
    matrix q;
    double u,alpha,w,t;
    m=a.mrow();
    n=a.mcol();
    if (m<n)
    return(0);
    q=eye(m);

    if(m==n)
    nn=m-1;
    else
    nn=n;

    for(k=0; k<=nn-1; k++)
    {
        u=0.0;
        for(i=k; i<=m-1; i++)
        {
            w=fabs(a(i,k));
            if(w>u)
            u=w;
        }
        alpha=0.0;
        for(i=k; i<=m-1; i++)
        {
            t=a(i,k)/u;
            alpha=alpha+t*t;
        }
        if( a(k,k)>0.0 )
        u=-u;
        alpha=u*sqrt(alpha);
        if( fabs(alpha)+1.0==1.0 )
        return(0);
        u=sqrt(2.0*alpha*(alpha-a(k,k)));
        if((u+1.0)!=1.0)
        {
            a(k,k)=(a(k,k)-alpha)/u;
            for(i=k+1; i<=m-1; i++)
            {
                a(i,k)=a(i,k)/u;
            }
            for(j=0; j<=m-1; j++)
            {
                t=0.0;
                for(jj=k; jj<=m-1; jj++)
                t=t+a(jj,k)*q(jj,j);
                for(i=k; i<=m-1; i++)
                //p=i*m+j;
                q(i,j)=q(i,j)-2.0*t*a(i,k);
              }
            for(j=k+1; j<=n-1; j++)
            {
                t=0.0;
                for(jj=k; jj<=m-1; jj++)
                t=t+a(jj,k)*a(jj,j);
                for(i=k; i<=m-1; i++)
                a(i,j)=a(i,j)-2.0*t*a(i,k);
              }
            a(k,k)=alpha;
            for(i=k+1; i<=m-1; i++)
            a(i,k)=0.0;
        }
    }
    return(1);
}
int qr1(matrix& a)
{
    int m=mrow(a);
    int n=mcol(a);
    int s;
    s=0;
    for(int i=0;i<n;i++)
    {
        vector w(m-s);
        vector u(m-s);
        double alpha;
        for(int j=s;j<m;j++)
        u[j-s]=a(j,i);
        double d=dot(u,u);
        alpha=sqrt(d);
        if(fabs(alpha)<0.00000001)
        {
            for(int j=s;j<m;j++)
            a(j,i)=0.0;
            continue;
        }
        if( a(s,i)>0.0 )
        alpha=-alpha;

        double r=sqrt(2.0*alpha*(alpha-u[0]));
        if(r>0.00000001)
        {
            u/=r;
            u[0]-=alpha/r;
            a(s,i)=alpha;
            for(j=s+1;j<m;j++)
            a(j,i)=0.0;
            for(j=i+1;j<n;j++)
            {
                d=0;
                for(int k=s;k<m;k++)
                d+=u[k-s]*a(k,j);
                for( k=s;k<m;k++)
                a(k,j)-=2.0*u[k-s]*d;
            }
        }
        s++;
    }
	return 1;
}

double rndnorm(double mu,double sigma)
{
    double m=(double)RAND_MAX;
    double data=0,n;
    for(int i=0;i<12;i++)
    {
        n=((double)rand())/m;
        data+=n;
    }
    n=(data-6)*sigma+mu;
    return n;
}
double rnduniform(double dmin, double dmax)
{
    double m=(double)RAND_MAX;
    double data=rand();
	data/=m;
    m=data*(dmax-dmin)+dmin;
    return m;
}
double cauchyrnd()
{
    double m=(double)RAND_MAX;
    double data=rand();
	data/=m;
    m=3.1415926*(data-0.5);
	data=tan(m);
    return data;
}


vector polyfit(vector& x,vector& y,int aicflag,int n)
{
    vector v,u;
    matrix a;
    int i,j;
    if(x()!=y())
    return v;
    double d;
    double delta,aicmin,newaic;
    if(aicflag<0)
    {
        a.initial(n,x());
        a[0]=1;
        for(i=1;i<n;i++)
        for(j=0;j<x();j++)
        a(i,j)=pow(x[j],i);
        SolveREqa(a,y,v);
    }
    else
    {
        if(aicflag)
        d=log(x())/x();
        else
        d=2.0/x();
        aicmin=10000000.0;
        for(int m=2;m<=n;m++)
        {
            a.initial(m,x());
            a[0]=1;
            for(i=1;i<m;i++)
            for(j=0;j<x();j++)
            a(i,j)=pow(x[j],i);
            delta=SolveREqa(a,y,u);
            delta=delta*delta/x();
            newaic=log(delta)+d*m;
            if(newaic<=aicmin)
            {
                aicmin=newaic;
                v=u;
            }
        }
    }
    return v;
}
vector polyval(vector& p,vector& x)
{
    int m=vlen(p);
    vector y;
    if( (x()==0) || (m==0))