#include<stdio.h>
#include<math.h>
intm //记录约束条件方程组的个数
intn //记录未知量的个数
floatM=1000000.0
floatA[100][100] //用于记录方程组的数目和系数
floatC[100] //用于存储目标函数中各个变量的系数
floatb[100] //用于存储常约束条件中的常数
floatCB[100]//用于存储基变量的系数
floatseta[100] //存放出基与入基的变化情况
floatcn[100] //存储检验数矩阵
floatx[100]
intnum[100] //用于存放出基与进基变量的情况
floatZ=0 //记录目标函数值
voidshuru()
voidprint()
intmincz()
intfind_line(int a)
voidexchange(int a,int b)
intmain()
{
int i,j=0
int p,q,temp //q：换入，拦含桐p：换出
shuru()
printf("\n--------------------------------------------------------------------------\n")
printf(" \tCB\tXB\tb\t")
for(i=0i<ni++)
printf(" X(%d)\t",i+1)
for(i=0i<ni++)
x[i]=0
printf("\n")
while(1) {
q=mincz()
if(q==-1) {
print()
printf("\n所得解已经是最优解!\n")
printf("\n最优解为:\n")
for(j=0j<mj++) {
temp=num[j]-1
x[temp]=b[j]
}
for(i=0i<ni++) {
printf("x%d=%.2f",i+1,x[i])
Z=Z+x[i]*C[i]
}
printf("Z=%.2f",Z)
break
}
print()
p=find_line(q)
printf("\np=%d,q=%d",p,q)
if(q==-1) break
exchange(p,q)
}
return 0
}
intmincz()
{
int i,k=0
int flag=0//检验数标记
float min=0
for(i=0i<ni++)
if(cn[i]>=0)
flag=1
else {
flag=0
break
}
if(flag==1)
return -1
//进行到此处，说明存在<0的检验数
//找到最小的检验数，作为换入变量
for(i=0i<ni++) {
if(min>cn[i]) {
min=cn[i]
k=i
}
}
return k
}
intfind_line(int a)
{
int i,k,j
int flag=0
float min
k=a
for(i=0i<mi++)
if(A[i][k]<=0)
flag=1
else {
flag=0
break
}
if(flag==1) {
printf("\n该线性规划无最优解!\n")
return -1
}
for(i=0i<mi++) {
if(A[i][k]>0)
seta[i]=b[i]/A[i][k]
else seta[i]=M
}
min=M
for(i=0i<mi++) {
if(min>=seta[i]) {
min=seta[i]
j=i
}
}
num[j]=k+1
CB[j]=C[k]
return j
}
voidexchange(int p,int q)
{
int i,j,c,l
float temp1,temp2,temp3
c=p//行号，换出
l=q//列号，换入
temp1=A[c][l] //A[c][l]主元
b[c]=b[c]/temp1
for(j=0j<nj++)
A[c][j]=A[c][j]/temp1 //主元化为1
for(i=0i<mi++) {
if(i!=c)
if(A[i][l]!=0) {
temp2=A[i][l]
b[i]=b[i]-b[c]*temp2
//主元所在列，其余元素化为0
for(j=0j<nj++)
A[i][j]=A[i][j]-A[c][j]*temp2
}
}
temp3=cn[l]
for(i=0i<ni++)
cn[i]=cn[i]-A[c][i]*temp3
}
voidprint()
{
int i,j=0
printf("\n--------------------------------------------------------------------------\n")
for(i=0i<mi++) {
printf("%8.2f\tX(%d) %8.2f",CB[i],num[i],b[i])
for(j=0j<nj++)
printf("%8.2f ",A[i][j])
printf("\n")
}
printf("\n--------------------------------------------------------------------------\n")
printf("\t\t\t")
for(i=0i<ni++)
printf(" %8.2f",cn[i])
printf("\n--------------------------------------------------------------------------\n")
}
voidshuru()
{
int i,j//循环变量
int k
printf("请输入线性规划问题的约束条件个数M：")
scanf("%d",&m)
printf("请输入线性规划问题的决策变量个数N：")
scanf("%d",&n)
printf("\n请输入方程组的系数矩阵A(%d行%d列):\n",m,n)
for(i=0i<mi++)
for(j=0j<nj++)
scanf("%f",&A[i][j])
printf("\n请输入初始基变量的数字代码矩阵:\n")
for(i=0i<mi++)
scanf("%d",&num[i])
printf("\n请输入方程组右边的值矩阵b:\n")
for(i=0i<mi++)
scanf("%f",&b[i])
printf("\n请输入目标函数各个变量的系数所构成的系数阵C:\n")
for(i=0i<ni++)
scanf("%f",&C[i])
for(i=0i<ni++)
cn[i]=-C[i]
for(i=0i<mi++) {
k=num[i]-1
CB[i]=C[k]
}
}
为方便查阅，再link一下 教材 。
假设我们讨论的LP问题有 个变量与 个限制条件。
基本概念部分大概就是这些。第一次看不用要求自己完全理解，不妨先继续往下学习再慢慢理解这些概念罩培。
单纯形方法的示意图如下：
我们先从一下四点直观的认识入手。
由于书上根本就没有优化 *** 作的概念，所以在此我自己给出其定义：
我们称增加的这个非基本变量为 输入变量 ，变为0的这个基本变量为 输出变量 。因为输入变量取代了输出变量成为新的基。（基的定义可以在前面找到。）
通俗一点地说就是选择一个非基本变量，让其不断增加，要求：
回到主题。之前提到单纯形法即对一个基本解实施若干次优化 *** 作后得到最优解的过程。我们先不考虑最优解的存在性，且断言： 任意非基本变量增加不使得目标函数增加等价于目标函数取得最大值。 这是因为由于符号的限制，非基本变量只能增加，而所有的变量都可以被非基本变量表示。因此非基本变量的增加包含了所有变量的各种变化。
当无法再进行优化 *** 作时有两种可能。一种是任意非基本变量的增加不使得态闷咐目标函数增加。此时目标函数取最大值。另外一种情况是任意非基本变量的增加不使得某一基本变量变为0。注意， LP问题的标准形式中对于变量的所有限制最终都会归结于符号限制 ，因此若基本变帆纯量不会变为0，则非基本变量可以无限地增加。 此时LP问题无界。
还有一点需要注意的是优化 *** 作一定会在有限步后结束，之后会讲到这一点。
经过上述讨论我们对单纯形算法的核心步骤应该有一个大致的了解了。
下面是一个小小的总结
单纯形算法到此结束
求多维函数极值的一种算法，由Nelder和Mead提出，又叫单纯形算法，但和线性规划中的单纯肢唤形算法是不同的，由于未利用任何求导运算，算法比较简单，但收敛速度较慢，适合变元数不是很多的方程求极值，算法的基本思想如下：给定n个特征，可以构造一个具有n+1个顶点的单纯形，初始化时需(n+1)*n维矩阵(看成是有n+1个顶点的单纯形) ，矩阵的每一行为n元向量，x0为第一行，xi=x0+r*ei，r为对问题的特征长度大小的估计值，ei为单位向量，x0可初始化为全为1的向量，即认为每个特征权重是相同的，然后选取其余的，在选取过程中，r可以取相同的值也可以取不同的值（r可以看作是对第i个特征权重的调整） 。
算法运行过程（以机器翻译中的rerank为例）：
假定BLEU=f(特征的和)，对n+1个顶点（n维向量）分别计算BLEU值（取相反数），然后从中选出BLEU(相反数)最大，次大和最小的三个点，算法每次都是把其中的最大点对应的各权重进行调整，使其变小向最小点靠拢，调整完毕后，计算其对应的BLEU，再从这些BLEU中选出BLEU(相反数)最大，次大和最小的三个点，一直迭代下去，直到最高点到最低点羡罩的比率范围合适或达到最大迭代次数为止。
源码：
double famoeb(double x[],vector<double>feat)
{//计算所有特征*权重的和
double y=0.0
for(int i=0i<FeatNumi++)
{
y+=x[i+1]*feat[i]
}
return y
}
//单纯形算法
void amoeba(double p[],double y[],int mp,int np,int ndim,double ftol,int&iter)
{
int i,j,ihi,inhi,mpts,nmax=20
double ypr,yprr,rtol,alpha=1.0
double beta=0.5
double gamma=2.0
int itmax=500
double pr[21],prr[21],pbar[21]
mpts=ndim+1
iter=0
do
{
int ilo=1
if(y[1]>y[2])
{
ihi=1
inhi=2
}
else
{
ihi=2
inhi=1
}
for(i=1i<=mptsi++)
{//寻找函数值中的最大，最小和次大值
if(y[i]<y[ilo])
{
ilo=i
}
if(y[i]>y[ihi])
{
inhi=ihi
ihi=i
}
else
{
if(y[i]>y[inhi])
{
if(i!=ihi)
{
inhi=i
}
}
}
}//结束寻找各种函数极值
rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]))//计算从最高点到最低点的比率范围，如合适则返回
if(rtol<ftol)
{
erase(pbar,prr,pr)
return
}
if(iter==itmax)//如到了最大的迭代次数，历派凯则返回
{
cout<<"amoeba exceeding maximum iterations."<<endl
return
}
iter=iter+1//进行下一次迭代
for(j=1j<=ndimj++)
{
pbar[j]=0.0
}
for(i=1i<=mptsi++)
{
if(i!=ihi)
{
for(j=1j<=ndimj++)
{
pbar[j]=pbar[j]+p[(i-1)*np+j]
}
}
}
for(j=1j<=ndimj++)
{
pbar[j]=pbar[j]/ndim
pr[j]=(1.0+alpha)*pbar[j]-alpha*p[(ihi-1)*np+j]//求反射点
}
vector<int>BestNo
ChooseOneBest(pr,numSentences,alldata,StartEndIndices,BestNo)
//开始计算BLEU值
vector<pairnum>initialScore(N_gram)
double referenceLength=0.0//参考翻译总长度
for(int k=0k<numSentencesk++)
{
int sent=BestNo[k]//当前句子的最好候选翻译的序号
for(int l=0l<N_graml++)
{
initialScore[l].left+=alldata[sent].ngram_data[l].left
initialScore[l].right+=alldata[sent].ngram_data[l].right
}
referenceLength+=alldata[sent].closest_length
}
ypr=-BLEU(initialScore,referenceLength)//计算本轮lamda所对应的bleu
if(ypr<=y[ilo])
{//得到一个比最佳点稍好的结果，用gamma做一次外推
for(j=1j<=ndimj++)
{
prr[j]=gamma*pr[j]+(1.0-gamma)*pbar[j]
}
vector<int>BestNo1
ChooseOneBest(prr,numSentences,alldata,StartEndIndices,BestNo1)
//开始计算BLEU值
vector<pairnum>initialScore1(N_gram)
double referenceLength1=0.0//参考翻译总长度
for(int m=0m<numSentencesm++)
{
int sent=BestNo1[m]//当前句子的最好候选翻译的序号
for(int n=0n<N_gramn++)
{
initialScore1[n].left+=alldata[sent].ngram_data[n].left
initialScore1[n].right+=alldata[sent].ngram_data[n].right
}
referenceLength1+=alldata[sent].closest_length
}
yprr=-BLEU(initialScore1,referenceLength1)//计算本轮lamda所对应的bleu
if(yprr<y[ilo])
{//以扩张点prr作为新的单纯形中的点
for(j=1j<=ndimj++)
{
p[(ihi-1)*np+j]=prr[j]
}
y[ihi]=yprr
}
else
{//以反射点pr作为单纯形中得点
for(j=1j<=ndimj++)
{
p[(ihi-1)*np+j]=pr[j]
}
y[ihi]=ypr
}
}
else
{//反射点不如最佳点，同次高点比较
if(ypr>=y[inhi])
{//反射点不如次高点，取一个中等程度低的点作一次一维收缩
if(ypr<y[ihi])
{
for(j=1j<=ndimj++)
{
p[(ihi-1)*np+j]=pr[j]
}
}
y[ihi]=ypr
for(j=1j<=ndimj++)
{
prr[j]=beta*p[(ihi-1)*np+j]+(1.0-beta)*pbar[j]
}
vector<int>BestNo2
ChooseOneBest(prr,numSentences,alldata,StartEndIndices,BestNo2)
//开始计算BLEU值
vector<pairnum>initialScore2(N_gram)
double referenceLength2=0.0//参考翻译总长度
for(int s=0s<numSentencess++)
{
int sent=BestNo2[s]//当前句子的最好候选翻译的序号
for(int t=0t<N_gramt++)
{
initialScore2[t].left+=alldata[sent].ngram_data[t].left
initialScore2[t].right+=alldata[sent].ngram_data[t].right
}
referenceLength2+=alldata[sent].closest_length
}
yprr=-BLEU(initialScore2,referenceLength2)//计算本轮lamda所对应的bleu
if(yprr<y[ihi])
{//以prr作为新单纯形中的点
for(j=1j<=ndimj++)
{
p[(ihi-1)*np+j]=prr[j]
}
y[ihi]=yprr//更新当前最高点出的函数值
}
else
{//单纯性太大，缩小原来的单纯形
for(i=1i<=mptsi++)
{
if(i!=ilo)
{
for(j=1j<=ndimj++)
{
pr[j]=0.5*(p[(i-1)*np+j]+p[(ilo-1)*np+j])
p[(i-1)*np+j]=pr[j]
}
vector<int>BestNo3
ChooseOneBest(pr,numSentences,alldata,StartEndIndices,BestNo3)
//开始计算BLEU值
vector<pairnum>initialScore3(N_gram)
double referenceLength3=0.0//参考翻译总长度
for(int u=0u<numSentencesu++)
{
int sent=BestNo3[u]//当前句子的最好候选翻译的序号
for(int v=0v<N_gramv++)
{
initialScore3[v].left+=alldata[sent].ngram_data[v].left
initialScore3[v].right+=alldata[sent].ngram_data[v].right
}
referenceLength3+=alldata[sent].closest_length
}
y[i]=-BLEU(initialScore3,referenceLength3)//计算本轮lamda所对应的bleu
}
}
}
}
else
{//反射点好于次高点，以反射点pr作为单纯形中得点
for(j=1j<=ndimj++)
{
p[(ihi-1)*np+j]=pr[j]
}
y[ihi]=ypr
}
}
}while(1)
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）