求助这个smo算法的

超级龟

这个程序中的拉格朗日乘子是怎么付值的，我怎么感觉都是0那？还有请高手给我以个训练文本，这个程序绝提实现的功能是什么？
#include "iostream.h"
#include "fstream.h"
#include "ctype.h"
#include "math.h"
#include "stdlib.h"
#include "time.h"

#define N 186//样本数总数
#define end_support_i 150
#define first_test_i 150
#define d 31//维数
//////////global variables
//int N=14104     //样本数
//int d=2    //维数
float C=100;//惩罚参数
float tolerance=0.001;//ui与边界0,C之间的公差范围
float eps=1e-3;//一个近似0的小数
float two_sigma_squared=10;//RBF核函数中的参数
float alph[end_support_i];//Lagrange multipiers
float b;//threshold
float error_cache[end_support_i];//存放non-bound样本误差
int target[N];//训练与测试样本的目标值
float precomputed_self_dot_product[N];//预存dot_product_func(i,i)的值，以减少计算量
int dense_points[N][d];//存放训练与测试样本，0-end_support_i-1训练;first_test_i-N测试

//函数的申明
int takeStep(int,int);
int examineNonBound(int);
int examineBound(int);
int examineFirstChoice(int,float);
int examineExample(int);
float kernel_func(int,int);
float learned_func(int);
float dot_product_func(int,int);
float error_rate();
void setX();
void setT();
void initialize();
///////////////////////////////////////////////////////
void main()
{
    ofstream os("d:\\smo1.txt");//存放训练的后的参数
    int numChanged=0,examineAll=1;
    //srand((unsigned int)time(NULL));
    initialize();
    //以下两层循环，开始时检查所有样本，选择不符合KKT条件的两个乘子进行优化，选择成功，返回1，否则，返回0
    //所以成功了，numChanged必然>0，从第二遍循环时，就不从整个样本中去寻找不符合KKT条件的两个乘子进行优化，
    //而是从边界的乘子中去寻找，因为非边界样本需要调整的可能性更大，边界样本往往不被调整而始终停留在边界上。
    //如果，没找到，再从整个样本中去找，直到整个样本中再也找不到需要改变的乘子为止，此时，算法结束。
    while(numChanged>0||examineAll)
    {
        numChanged=0;
        if(examineAll)
        {
            for(int k=0;k<end_support_i;k++)
                numChanged+=examineExample(k);//examin all exmples
        }
        else{
            for(int k=0;k<end_support_i;k++)
                if(alph[k]!=0&&alph[k]!=C)
                    numChanged+=examineExample(k);//loop k over all non-bound lagrange multipliers
            }
        if(examineAll==1)
            examineAll=0;
        else if(numChanged==0)
            examineAll=1;
    }
    //存放训练后的参数
    {
        os<<"d="<<d<<endl;
        os<<"b="<<b<<endl;
        os<<"two_sigma_squared="<<two_sigma_squared<<endl;
        int n_support_vectors=0;
        for(int i=0;i<end_support_i;i++)
            if(alph>0)//此地方是否可以改为alph>tolerance?????????????????
                n_support_vectors++;
        os<<"n_support_vectors="<<n_support_vectors<<"rate="<<(float)n_support_vectors/first_test_i<<endl;
        for(i=0;i<end_support_i;i++)
            if(alph>0)
            os<<"alph["<<i<<"]="<<alph<<" ";
        os<<endl;
    }
    //输出，以及测试
    cout<<error_rate()<<endl;
}

/////////examineExample程序
//假定第一个乘子ai(位置为i1)，examineExample(i1)首先检查，如果它超出tolerance而违背KKT条件，那么它就成为第一个乘子；
//然后，寻找第二个乘子(位置为i2),通过调用takeStep(i1,i2)来优化这两个乘子。
int examineExample(int i1)
{
    float y1,alph1,E1,r1;
    y1=target[i1];
    alph1=alph[i1];
    if(alph1>0&&alph1<C)
        E1=error_cache[i1];
    else
        E1=learned_func(i1)-y1;//learned_func为计算输出函数
    r1=y1*E1;
    //////违反KKT条件的判断
    if((r1>tolerance&&alph1>0)||(r1<-tolerance&&alph1<C))
    {
        /////////////使用三种方法选择第二个乘子
        //1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
        //2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
        //3：如果上面也失败，则从随机位置查找整个样本,改为bound样本
        if (examineFirstChoice(i1,E1))//1
             {
                     return 1;
              }

              if (examineNonBound(i1))//2
              {
                 return 1;
              }

              if (examineBound(i1))//3
             {
                 return 1;
              }
    }
    ///没有进展
          return 0;
}

////////////1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
int examineFirstChoice(int i1,float E1)
{
   int k,i2;
   float tmax;
   for(i2=-1,tmax=0.0,k=0;k<end_support_i;k++)//*******************************end_support_i
       if(alph[k]>0&&alph[k]<C)
       {
           float E2,temp;
           E2=error_cache[k];
           temp=fabs(E1-E2);
           if(temp>tmax)
           {
               tmax=temp;
               i2=k;
           }
       }
   if(i2>=0)
   {
       if(takeStep(i1,i2))
           return 1;
   }
   return 0;
}
///////////////    2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
int examineNonBound(int i1)
{
       int k,k0 = rand()%end_support_i;
        int i2;
       for (k = 0; k < end_support_i; k++)
       {
              i2 = (k + k0) % end_support_i;//从随机位开始

              if (alph[i2] > 0.0 && alph[i2] < C)
              {
                 if (takeStep(i1, i2))
                 {
                    return 1;
              }
              }
       }
       return 0;
}
////////////3：如果上面也失败，则从随机位置查找整个样本,(改为bound样本)
int examineBound(int i1)
{
       int k,k0 = rand()%end_support_i;
        int i2;
       for (k = 0; k < end_support_i; k++)
       {
              i2 = (k + k0) % end_support_i;//从随机位开始

              //if (alph[i2]= 0.0 || alph[i2]=C)//修改******************
              {
                 if (takeStep(i1, i2))
                 {
                    return 1;
              }
              }
       }
       return 0;   
}
///////////takeStep()
//用于优化两个乘子，成功，返回1，否则，返回0
int takeStep(int i1,int i2)
{
    int y1,y2,s;
    float alph1,alph2;//两个乘子的旧值
    float a1,a2;//两个乘子的新值
    float E1,E2,L,H,k11,k22,k12,eta,Lobj,Hobj,delta_b;
   
    if(i1==i2) return 0;//不会优化两个同一样本
    //给变量赋值
    alph1=alph[i1];
    alph2=alph[i2];
    y1=target[i1];
    y2=target[i2];
    if(alph1>0&&alph1<C)
        E1=error_cache[i1];
    else
        E1=learned_func(i1)-y1;//learned_func(int)为非线性的评价函数，即输出函数
    if(alph2>0&&alph2<C)
        E2=error_cache[i2];
    else
        E2=learned_func(i2)-y2;
    s=y1*y2;
    //计算乘子的上下限
    if(y1==y2)
    {
        float gamma=alph1+alph2;
        if(gamma>C)
        {
            L=gamma-C;
            H=C;
        }
        else
        {
            L=0;
            H=gamma;
        }
    }
    else
    {
        float gamma=alph1-alph2;
        if(gamma>0)
        {
            L=0;
            H=C-gamma;
        }
        else
        {
            L=-gamma;
            H=C;
        }
    }//计算乘子的上下限
    if(L==H) return 0;
    //计算eta
    k11=kernel_func(i1,i1);//kernel_func(int,int)为核函数
    k22=kernel_func(i2,i2);
    k12=kernel_func(i1,i2);
    eta=2*k12-k11-k22;
    if(eta<0)
    {
        a2=alph2+y2*(E2-E1)/eta;//计算新的alph2
        //调整a2，使其处于可行域
        if(a2<L)     a2=L;
        if(a2>H)    a2=H;
    }
    else//此时，得分别从端点H,L求目标函数值Lobj,Hobj，然后设a2为求得最大目标函数值的端点值
    {
        float c1=eta/2;
        float c2=y2*(E2-E1)-eta*alph2;
        Lobj=c1*L*L+c2*L;
        Hobj=c1*H*H+c2*H;
        if(Lobj>Hobj+eps)//eps****************
            a2=L;
        else if(Lobj<Hobj-eps)
            a2=H;
        else
            a2=alph2;//加eps的目的在于，使得Lobj与Hobj尽量分开，如果，很接近，就认为没有改进(make progress)
    }
    if(fabs(a2-alph2)<eps*(a2+alph2+eps))
        return 0;
    a1=alph1-s*(a2-alph2);//计算新的a1
    if(a1<0)//调整a1,使其符合条件*****??????????????????????????????????????????
    {
        a2+=s*a1;
        a1=0;
    }
    else if(a1>C)
    {
        float t=a1-C;
        a2+=s*t;
        a1=C;
    }
    //更新阀值b
    {
        float b1,b2,bnew;
        if(a1>0&&a1<C)
            bnew=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
        else{
            if(a2>0&&a2<C)
                bnew=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
            else{
                b1=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
                b2=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
                bnew=(b1+b2)/2;
            }
        }
        delta_b=bnew-b;
        b=bnew;
    }
   
    //对于线性情况，要更新权向量，这里不用了
    //更新error_cache，对取得进展的a1,a2,所对应的i1,i2的error_cache[i1]=error_cache[i2]=0
    {
        float t1=y1*(a1-alph1);
        float t2=y2*(a2-alph2);
        for(int i=0;i<end_support_i;i++)
        if(0<alph&&alph<C)
            error_cache+=t1*kernel_func(i1,i)+t2*(kernel_func(i2,i))-delta_b;
            error_cache[i1]=0.0;   
            error_cache[i2]=0.0;
    }
    alph[i1]=a1;//store a1,a2 in the alpha array
    alph[i2]=a2;
    return 1;//说明已经取得进展
}

//////////learned_func(int)
//评价分类学习函数
float learned_func(int k)
{
    float s=0.0;
    for(int i=0;i<end_support_i;i++)
        if(alph>0)
            s+=alph*target*kernel_func(i,k);
        s-=b;
    return s;
}
//计算点积函数dot_product_func(int,int)
float dot_product_func(int i1,int i2)
{
    float dot=0;
    for(int i=0;i<d;i++)
        dot+=dense_points[i1]*dense_points[i2];   
    return dot;
}
//径向机核函数RBF:kernel_func(int,int)
float kernel_func(int i1,int i2)
{
    float s=dot_product_func(i1,i2);
    s*=-2;
    s+=precomputed_self_dot_product[i1]+precomputed_self_dot_product[i2];  
    return exp(-s/two_sigma_squared);
}
//初始化initialize()
void initialize()
{
    //初始化阀值b为0
    b=0.0;
    //初始化alph[]为0
    for(int i=0;i<end_support_i;i++)
    alph=0.0;
   
    //设置样本值矩阵
    setX();
    //设置目标值向量
    setT();
    //设置预计算点积
    {for(int i=0;i<N;i++)
        precomputed_self_dot_product=dot_product_func(i,i);//注意，这是对训练样本的设置，对于测试样本也应考虑?????????????????
    }
}
//计算误差率error_rate()
float error_rate()
{    ofstream to("d:\\smo_test.txt");
    int tp=0,tn=0,fp=0,fn=0;
    float ming=0,te=0,total_q=0,temp=0;
    for(int i=first_test_i;i<N;i++)
    {    temp=learned_func(i);
        if(temp>0&&target>0)
            tp++;
        else if(temp>0&&target<0)
            fp++;
        else if(temp<0&&target>0)
            fn++;
        else if(temp<0&&target<0)
            tn++;
        to<<i<<"  实际输出"<<temp<<endl;
    }
    total_q=(float)(tp+tn)/(float)(tp+tn+fp+fn);//总精度
    ming=(float)tp/(float)(tp+fn);
    te=(float)tp/(float)(tp+fp);
    to<<"---------------测试结果-----------------"<<endl;
    to<<"tp="<<tp<<"   tn="<<tn<<"  fp="<<fp<<"  fn="<<fn<<endl;
    to<<"ming="<<ming<<"  te="<<te<<"  total_q="<<total_q<<endl;
    return (1-total_q);
}
//计算样本X[]
void setX(){
    ifstream ff("D:\\17_smo.txt");////////////////////////
    if(!ff)
    {
        cerr<<"error!!"<<endl;
        exit(1);
    }
    int i=0,j=0;
    char ch;
    while(ff.get(ch))
    {
        if(isspace(ch)) continue;
        if(isdigit(ch))
        {
            dense_points[j]=(int)ch-48;
            j++;
            if(j==d)
            {
                j=0;
                i++;
            }
        }
    }
}
//set targetT[]
void setT()
{
    for(int i=0;i<75;i++)
        target=1;
    for(i=75;i<150;i++)
        target=-1;
    for(i=150;i<168;i++)
        target=1;
    for(i=168;i<186;i++)
        target=-1;
}

荒地有群狼

我也在烦这个程序呢，楼主研究出来的话，告诉声