ant.txt

C++的蚂蚁算法源代码

  for(i=0;   i<j;   i++)   
  a[i+j*n]=0.0;   
  return   md;   
  }   
    
  DOUBLE   lineopt(matrix&   aa,matrix&   bb,   matrix&   cc,   matrix   &   xx)   
    //   线性规划最优点寻找程序，aa为mXn不等式约束条件左端系数矩阵，bb为不等式约束   
    //   条件的右端项，为m维向量，cc为目标函数系数，n维向量，xx返回极小点，为n维向量   
  {   
  if(aa.rownum   !=   bb.rownum   ||   aa.colnum   !=   cc.rownum   ||   
  aa.colnum   !=   xx.rownum)   throw   TMESSAGE("dimenstion   not   right!");   
  size_t   n=aa.colnum,   m=aa.rownum;   
  size_t   i,mn,k,j;   
  matrix   a(m,n+m);   
  for(i=0;i<m;i++)   {   
  for(j=0;j<n;j++)   
  a.set(i,j,aa(i,j));   
  for(j=n;j<n+m;j++)   
  if(j-n   ==   i)   a.set(i,j,1.0);   
  else   a.set(i,j,0.0);   
  }   
  matrix   c(m+n);   
  c   =   0.0;   
  for(i=0;i<m;i++)   
  c.set(i,cc(i));   
  lbuffer*   jjs   =   getnewlbuffer(m);   
  lbuffer&   js   =   (*jjs);   
  DOUBLE   s,z,dd,y;   //,*p,*d;   
    
  for   (i=0;   i<m;   i++)   js[i]=n+i;   
  matrix   p(m,m);   
  matrix   d;   
  mn=m+n;   
  matrix   x(mn);   
  while   (1)   
  {   for   (i=0;   i<m;   i++)   
    for   (j=0;   j<m;   j++)   
  p.set(i,j,a(i,(size_t)js[j]));   
      p.inv();   
  d   =   p*a;   
  x   =   0.0;   
      for   (i=0;   i<m;   i++)   
    {   s=0.0;   
  for   (j=0;   j<=m-1;   j++)   
  s+=p(i,j)*bb(j);   
  x.set((size_t)js[i],s);   
    }   
      k=mn;   dd=1.0e-35;   
      for   (j=0;   j<mn;   j++)   
    {   z=0.0;   
  for   (i=0;   i<m;   i++)   
  z+=c((size_t)js[i])*d(i,j);   
  z-=c(j);   
  if   (z>dd)   {   dd=z;   k=j;}   
    }   
      if   (k==mn)   
    {   s=0.0;   
  for   (j=0;   j<n;   j++)   {   
  xx.set(j,x(j));   
  s+=c(j)*x(j);   
  }   
  delete   jjs;   
  return   s;   
    }   
      j=m;   
      dd=1.0e+20;   
      for   (i=0;   i<=m-1;   i++)   
    if   (d(i,k)>=1.0e-20)       //   [i*mn+k]>=1.0e-20)   
  {   y=x(size_t(js[i]))/d(i,k);   
      if   (y<dd)   {   dd=y;   j=i;}   
  }   
      if   (j==m)   {   delete   jjs;   
  throw   TMESSAGE("lineopt   failed!");   
  }   
      js[j]=k;   
  }   
  }   
  matrix   matrix::mean(void){   
      matrix   mat(1,getcolnum());   
      double   sum,rownum;   
      rownum=(double)getrownum();   
      for(int   j=0;j<(int)getcolnum();j++){   
          sum=0.0;   
          for(int   i=0;i<(int)getrownum();i++)   
              sum+=value(i,j);   
          mat.set(0,j,sum/rownum);   
      }   
      return   mat;   
  }   
  matrix   matrix::cov(void){   
      matrix   mat_cov(getcolnum(),getcolnum());   
      matrix   mat_mean;   
      mat_mean=mean();   
      double   sum,rownum;   
      rownum=(double)getrownum();   
      for(int   i=0;i<(int)getcolnum();i++)   
          for(int   j=i;j<(int)getcolnum();j++){   
              sum=0.0;   
              for(int   k=0;k<(int)getrownum();k++)   
                  sum+=value(k,i)*value(k,j);   
              sum-=(rownum*mat_mean(0,i)*mat_mean(0,j));   
              sum/=(rownum-1);   
              mat_cov.set(i,j,sum);   
              mat_cov.set(j,i,sum);   
          }   
      return   mat_cov;   
  }   
  #endif   //MATRIX_CPP   
    
    
  文件4：matrix.h   
    
  #ifndef   MATRIX_H   
    
  #define   MATRIX_H   
    
  //   定义适合数学运算的实数matrix类，数据将存放在buffer类中   
    
  #include   <iostream.h>   
  #include   <iomanip.h>   
  #include   <stdio.h>   
    
  //   buffer.h包含实数缓存类buffer的定义   
  #include   "buffer.cpp"   
    
  //   实数矩阵类，每个矩阵有rownum行，colnum列，总共rownum乘colnum个实数   
  //   下标一律从0开始起算，即存在第0行第0列   
  class   matrix   {   
    public:   
  buffer   *   buf;     //   缓存区指针，指向存放实数的地方，可以是任何媒体   
  //   buffer类是一抽象类，具体的子类可以将数据存放在内存或者   
  //   临时文件中，用户还可以设计自己的缓存区子类   
  size_t   rownum,   colnum;     //   矩阵的行数和列数   
  unsigned istrans:1; //   转置标志，0为缺省，1为转置   
  unsigned   isneg:1; //   为负标志，0为非负，1为取负   
  unsigned   issym:1; //   对称标志，0为非对称，1为对称，不具有约束力   
  //   在下面的所有构造函数中，缓存区指针b通常不设，缺省为0，这样程序将根据   
  //   自动申请缺省的缓存区类型。如果用户定义了特殊的缓顾虑区类型，可以   
  //   先申请后在构造函数中指明   
  matrix(buffer   *   b=0);   //   缺省构造函数，产生0行0列空矩阵   
  matrix(size_t   n,   buffer   *   b=0);   //   产生nX1列向量   
  matrix(size_t   r,   size_t   c,   buffer   *   b=0);   //   构造函数，产生r行c列矩阵   
  matrix(matrix&   m);   //拷贝构造函数，产生的矩阵内容同m相同，   
  //   新产生的矩阵并不申请新的缓存，而是指向与m同样的缓存，   
  //   只是它们共同指向的缓存的引用数增1，这样可以节省存储资源   
  matrix(const   char   *   filename,   buffer   *   b=0);   /*   文件构造函数，   
  使用数据文件对矩阵初使化，数据文件为文本文件，格式：   
  以小写的matrix关键词打头，然后是空格，再跟矩阵的行数，再跟   
  空格，再跟列数，再跟空格，接着是内容，   
  内容按行分，首先是行号，从1开始计数，行号后面紧跟着冒号：，   
  然后是这行的数据，各个数据之间都用空格或者换行隔开   */   
  matrix(void   *   data,   size_t   r,   size_t   c=1,   buffer   *   b=0);   /*   数据构造函数，   
  行数为r,列数为c,data必须指向一存放实数的内存，按先列后行的   
  顺序进行存放。构造函数进行数据拷贝，因此在构造函数完成后,   
  data指向的数据区可以释放   */   
        size_t   getrownum(){return   istrans?colnum:rownum;   };//获得矩阵的行数   
        size_t   getcolnum(){return   istrans?rownum:colnum;   };//获得矩阵的列数   
        void   loadfile(const   char   *   filename);//从指定数据文件中获得矩阵数据   
        void   savefile(char   *   filename,int   g=2);//存矩阵数据到指定数据文件中   
  DOUBLE   &   operator()(size_t   r,   size_t   c);   //   重载函数运算符()返回第r行第c列的实数，   
  DOUBLE   &   operator()(size_t   r){return   (*this)(r,0);};   //   取回第r行0列的实数   
  virtual   void   set(size_t   r,   size_t   c,   DOUBLE   v)   {   
  size_t   l;   
  v   =   isneg   ?   -v   :   v;   
  l   =   istrans   ?   r+c*rownum   :   r*colnum+c;   
  buf   =   buf->set(l,   v);   };   //   设置第r行第c列的值为v，   
  //   在设置时如果缓存的引用数大于1，产生新的缓存   
  //   因此每次设置返回的缓存区指针都要再赋给buf。   
  void   set(size_t   r,   DOUBLE   v)   {   
  set(r,0,v);}; //   设置第r行0列的值为v,用作列向量   
  virtual   ~matrix(){   //   析构函数，如果缓存的引用数大于1，说明还有其它的   
  //   矩阵在用此缓存，则不释放缓存，只是将引用数减1   
  buf->refnum--;   
  if(!buf->refnum)   delete   buf;};   
  matrix   operator=(matrix   m);   //   重载赋值运算符，使矩阵的内容等于另一个矩阵   
  matrix   operator=(DOUBLE   a);   //   通过赋值运算符将矩阵所有元素设为a   
        matrix   operator-();   
  matrix   operator*(DOUBLE   a); //   重载乘法运算符，实现矩阵与数相乘，矩阵在前   
  //   数在后   
  matrix   operator*=(DOUBLE   a);     //   重载*=运算符，功能同上，但*产生了一个新   
  //   的矩阵，原矩阵不变，而本运算符修改了原矩阵   
  friend   matrix   operator*(DOUBLE   a,   matrix&   m);   
  //   数乘矩阵，但数在前矩阵在后   
  matrix   operator+(matrix&   m);   //   矩阵相加，原二矩阵不变并产生一个新的和矩阵   
  matrix   operator+=(matrix   m);   //   矩阵相加，不产生新矩阵，有一矩阵的内容改变为和   
    
  matrix   operator+(DOUBLE   a); //   矩阵加常数，指每一元素加一固定的常数，产生   
  //   新矩阵，原矩阵不变   
  friend   matrix   operator+(DOUBLE   a,   matrix&   m);   //   常数加矩阵，产生新的和矩阵   
  matrix   operator+=(DOUBLE   a); //   矩阵自身加常数，自身内容改变   
    
  matrix   operator*(matrix   m);   //   矩阵相乘，产生出新的矩阵   
  matrix   operator*=(matrix   m);   //   矩阵相乘，结果修改原矩阵   
    
  //matrix   operator-();   //   矩阵求负，产生新的负矩阵   
  matrix     trans(){   
  size_t   r   =   rownum;   rownum   =   colnum;   colnum   =   r;   
  istrans   =   !istrans;   return   (*this);};   //   矩阵自身转置   
  matrix   t();   //   矩阵转置，产生新的矩阵   
  matrix     neg(){isneg   =   !isneg;   return   (*this);};   //   将自己求负   
    
  matrix   operator-(matrix   m);   //   矩阵相减，产生新的矩阵   
  matrix   operator-=(matrix   m);   //   矩阵相减，结果修改原矩阵   
  matrix   operator-(DOUBLE   a)   //   矩阵减常数，指每一元素减一固定的常数，产生   
  //   新矩阵，原矩阵不变   
  {   return   (*this)+(-a);   };   
  friend   matrix   operator-(DOUBLE   a,   matrix   m);   //   常数减矩阵，产生新的矩阵   
  matrix   operator-=(DOUBLE   a) //   矩阵自身减常数，自身内容改变   
  {   return   operator+=(-a);};   
    
  int   isnear(matrix   m,   double   e   =   defaulterr);   //   检查两个矩阵是否近似相等,   
  //   如近似相等则返回1，否则返回0，e为允许误差   
  int   isnearunit(double   e   =   defaulterr);   //   检查矩阵是否近似为单位矩阵   
  //   如是则返回1，否则返回0，e为允许误差   
    
  matrix   row(size_t   r);   //   提取第r行行向量   
  matrix   col(size_t   c);   //   提取第c列列向量   
    
        matrix   mean(void);   
        matrix   cov(void);   
    
  void   checksym(); //   检查和尝试调整矩阵到对称矩阵   
    
  void   swapr(size_t   r1,   size_t   r2,   size_t   k=0);   //   交换两行,只交换k列以上   
  void   swapc(size_t   c1,   size_t   c2,   size_t   k=0);   //   交换两列,只交换k行以上   
  void   swap(size_t   r1,   size_t   c1,   size_t   r2,   size_t   c2){   //   交换两个元素   
  DOUBLE   a;   a=(*this)(r1,c1);set(r1,c1,(*this)(r2,c2));   
  set(r2,c2,a);};   
  DOUBLE   maxabs(size_t   &r,   size_t   &c,   size_t   k=0);   
    //   求主元值和位置，即k行k列之后的最大元素，   
  //   元素所在的行列将放在r,c中，返回此元素值   
  size_t   zgsxy(matrix   &   m,   int   fn=0); //   主高斯消元运算   
  matrix   operator/=(matrix   m);   //   用主高斯消元法求解线性方程的解   
  //   矩阵本身为系数矩阵，m为常数向量，返回解向量   
  matrix   operator/(matrix   m);   //   左乘m的逆矩阵   
    
  matrix   inv(); //   用全选主元高斯-约当法求逆矩阵   
  matrix   operator~(); //   求逆矩阵，但产生新矩阵   
  friend   matrix   operator/(DOUBLE   a,   matrix   m);   //   求逆矩阵再乘常数   
  matrix   operator/=(DOUBLE   a);   //   所有元素乘a的倒数，自身改变   
  matrix   operator/(DOUBLE   a);   //   所有元素乘a的倒数，产生新的矩阵   
  matrix   cholesky(matrix   &d); //   用乔里斯基分解法求对称正定阵的线性   
  //   方程组ax=d返回方程组的解，本身为a，改变为分解阵u,d不变   
    
  DOUBLE   det(DOUBLE   err=defaulterr); //   求行列式的值   
  size_t   rank(DOUBLE   err=defaulterr); //   求矩阵的秩   
    
  void   house(buffer   &   b,   buffer   &   c);   //   用豪斯荷尔德变换将对称阵变为对称三对   
  //   角阵，b返回主对角线元素，c返回次对角线元素   
  void   trieigen(buffer&   b,   buffer&   c,   size_t   l=600,   DOUBLE   eps=defaulterr);   
  //   计算三对角阵的全部特征值与特征向量   
  matrix   eigen(matrix   &   eigenvalue,   DOUBLE   eps=defaulterr,   size_t   l=600);   
    //   计算对称阵的全部特征向量和特征值   
  //   返回按列排放的特征向量，而eigenvalue将返回一维矩阵，为所有特征值   
  //   组成的列向量   
  int   ispositive(); //   判定矩阵是否对称非负定阵，如是返回1，否则返回0   
    
  void   hessenberg(); //   将一般实矩阵约化为赫申伯格矩阵   
  void   qreigen(matrix   &   a,   matrix   &   b,   size_t   l=600,   DOUBLE   eps=defaulterr);   
  //   求一般实矩阵的所有特征根   
  //   a和b均返回rownum行一列的列向量矩阵，返回所有特征根的实部和虚部   
    
  friend   ostream&   operator<<(ostream&   o,   matrix&   m);   
  friend   istream&   operator>>(istream&   in,   matrix&   m);   
        DOUBLE   &   value(size_t   r,   size_t   c){   //   返回第r行c列的数值，不做行列检查   
            DOUBLE   &a   =   istrans   ?   (*buf)[r+c*rownum]   :   (*buf)[r*colnum   +   c];   
  a=isneg   ?   -a   :   a;   
              return   a;   
  };   
    protected:   
  DOUBLE   detandrank(size_t   &   r,   DOUBLE   err); //   求方阵的行列式和秩   
  };   
    
  inline   matrix   operator*(DOUBLE   a,   matrix&   m)   {   
  return   m*a;   
  };   
    
  inline   matrix   operator+(DOUBLE   a,   matrix&   m)   {   //   常数加矩阵，产生新的和矩阵   
  return   m+a;   
  };   
    
  matrix   operator-(DOUBLE   a,   matrix&   m);   //   常数减矩阵，产生新的矩阵   
  matrix   operator/(DOUBLE   a,   matrix&   m);   //   求逆矩阵再乘常数   
    
    
  ostream&   operator<<(ostream&   o,   matrix&   m);   
  istream&   operator>>(istream&   in,   matrix&   m);   
    
  matrix   unit(size_t   n); //   产生n阶单位矩阵   
    
  DOUBLE   gassrand(int   rr=0); //   返回一零均值单位方差的正态分布随机数   
  //   rr为种子   
    
  class   gassvector   //返回零均值给定协方差阵的正态随机向量的类   
  {   
    public:   
  matrix   a; //   a是增益矩阵，由协方差矩阵算出   
  gassvector(matrix   &   r); //r必须是正定对称阵，为正态随机向量的协方差   
  matrix   operator()(matrix   &   r);   //   返回给定协方差矩阵的正态随机向量   
  matrix   operator()(); //   返回已设定的协方差矩阵的正态随机向量   
  void   setdata(matrix   &   r);   //   根据协方差矩阵设置增益矩阵   
  };   
    
  DOUBLE   lineopt(matrix&   a,matrix&   b,   matrix&   c,   matrix   &   x);   
  //   线性规划最优点寻找程序，a为mXn不等式约束条件左端系数矩阵，b为不等式约束   
  //   条件的右端项，为m维向量，c为目标函数系数，n维向量，x返回极小点，为n维向量   
    
  #endif   //   MATRIX_H   
    
  文件5：buffer.cpp   
    
  #ifndef   BUFFER_CPP   
  #define   BUFFER_CPP   
  #include   <stdio.h>   
  #include   <string.h>   
  #include   <fstream.h>   
  #include   <math.h>   
    
  //   功能：实现buffer.h中的各个缓冲类的函数   
    
  #include   "buffer.h"   
    
  DOUBLE   defaulterr   =   1e-8;   //   缺省的误差调整值   
  int   doadjust   =   1; //   决定数据是否根据误差值自动向整数靠扰，   
  //   其实是向小数点后三位靠扰   
    
  void   membuffer::Set(size_t   n,   DOUBLE   d)   
  {   
  //   如调整标志设置，则调整   
  if(doadjust)   adjust(d);   
  if(refnum   ==   1)   { //   如果引用数为1，则将第n个实数的值设为d，并返回自己   
  //   的指针   
  retrieve(n)   =   d;   
  return   this;   
  }   
  refnum   --; //   引用数大于1，因此将本缓冲区的引用数减1   
  buffer   *   b;   
  b   =   clone(); //   克隆一个内容同自己一样的新缓冲区   
  (b->retrieve(n))   =   d;   //