#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/timeb.h>
#define PI1 0
#define PI2 0
#define PI3 3
#define PI4 0
void intswap(int *a, int *b)
{
int temp;
temp = *a;
*a = *b;
*b = temp;
}
void swap(double *a, double *b)
{
double temp;
temp = *a;
*a = *b;
*b = temp;
}
void *MySort(double *list, int intNumofVars,int *Index)
{
int i,j,iMax;
for(i=0;i<intNumofVars;i++)
{
iMax=i;
for(j=i+1;j<intNumofVars;j++)
{
if(list[j]>list[iMax])
{
iMax=j;
}
}
if(iMax!=i)
{
swap(&list[i],&list[iMax]);
intswap(&Index[i],&Index[iMax]);
}
}
return Index;
}
double *ContinuousKP(int *Pi01, int *W, int c, int n, double *B)
{
double *Xcont=NULL;
int *Xbin=NULL;
int *Index=NULL;
double *Pi0=NULL;
int ind=0;
int resC;
double UB1=0;
double LB1=0;
double ma=0;
int i,k;
resC=c;
Xcont=(double*)calloc(n,sizeof(double));
Xbin=(int*)calloc(n,sizeof(int));
Index=(int*)calloc(n,sizeof(int));
Pi0=(double*)calloc(n,sizeof(double));
for(i=0;i<n;i++)
{
Pi0[i]=(double) Pi01[i]/W[i];
}
for(i=0;i<n;i++)
{
ma=Pi0[0];
ind=0;
for(k=0;k<n;k++)
{
if(Pi0[k]>ma)
{
ma=Pi0[k];
ind=k;
}
}
Index[i]=ind;
if(resC>=W[ind])
{
Xcont[ind]=1;
Xbin[ind]=1;
UB1=UB1+Pi01[ind];
LB1=LB1+Pi01[ind];
resC=resC-W[ind];
Pi0[ind]=-32000;
}
}
if(resC>0)
{
ma=Pi0[0];
ind=0;
for(k=0;k<n;k++)
{
if(Pi0[k]>ma)
{
ma=Pi0[k];
ind=k;
}
}
Xcont[ind]=(double) resC/W[ind];
UB1=UB1+ Xcont[ind]*Pi01[ind];
}
B[0]=UB1;
B[1]=LB1;
free(Xcont);
free(Xbin);
free(Index);
free(Pi0);
return B;
}
int *FillUpExchange(int *y, int **Profit, int *W, int c, int n)
{
int i,j,k,l;
int cy=0, cdum=0;
int *dummy=NULL;
int valdum=0;
dummy = (int*)calloc(n,sizeof(int));
/* Fill up*/
for(i=0;i<n;i++)
{
cy+=y[i]*W[i];
}
for(i=0;i<n;i++)
{
if(y[i] == 0 && cy + W[i] <= c)
{
y[i] = 1;
cy += W[i];
}
}
/* Exchange*/
for(i=0;i<n-1;i++)
{
for(j=i+1;j<n;j++)
{
if(y[i] == 1 && y[j] == 0)
{
for(k=0;k<n;k++)
{
dummy[k] = y[k];
}
dummy[i] = 0;
dummy[j] = 1;
for(k=0;k<n;k++)
{
cdum = cdum + dummy[k]*W[k];
for(l=0;l<n;l++)
{
valdum += y[k]*Profit[k][l]*y[l] - dummy[k]*Profit[k][l]*dummy[l];
}
}
if(cdum <= c && valdum <= 0)
{
y[i] = 0;
y[j] = 1;
}
}
else if(y[i] == 0 && y[j] == 1)
{
for(k=0;k<n;k++)
{
dummy[k] = y[k];
}
dummy[i] = 1;
dummy[j] = 0;
for(k=0;k<n;k++)
{
cdum = cdum + dummy[k]*W[k];
for(l=0;l<n;l++)
{
valdum += y[k]*Profit[k][l]*y[l] - dummy[k]*Profit[k][l]*dummy[l];
}
}
if(cdum <= c && valdum <= 0)
{
y[i] = 1;
y[j] = 0;
}
}
valdum = 0;
cdum = 0;
}
}
free(dummy);
return y;
}
int *Dynamic(int *p,int *w,int c, int n,int *x)
{
/* This routine solves the linear KP of the form
max px
s.t wx<=c
x_i \in {0,1} for i=1,...,n.
The profit vector p, the weight vector w, and the
capacity c are made of integer numbers. And n is the
number of items.
p and w should be augmented to have n+1 elements.
So p=[0,p] and w=[0,w].
*/
int k,r,i,j;
int ma, rstar;
int **f=NULL;
int *fend=NULL;
f=(int**)calloc(n+1,sizeof(int*));
fend=(int*)calloc(c+1,sizeof(int));
for(i=0;i<n+1;i++)
{
f[i]=(int*)calloc(c+1,sizeof(int));
for(j=0;j<c+1;j++)
{
f[i][j]=-30000;
}
}
f[0][0]=0;
for(k=0;k<n;k++)
{
for(r=0;r<c+1;r++)
{
if(f[k][r]>f[k+1][r])
{
f[k+1][r]=f[k][r];
}
if(r+w[k+1]<=c && f[k][r]+p[k+1]>f[k+1][r+w[k+1]])
{
f[k+1][r+w[k+1]]=f[k][r]+p[k+1];
}
}
}
for(i=0;i<c+1;i++)
{
fend[i]=f[n][i];
}
ma=fend[0];
rstar=0;
for(i=0;i<c+1;i++)
{
if(fend[i]>ma)
{
ma=fend[i];
rstar=i;
}
}
for(k=n;k>0;k--)
{
if(f[k][rstar]!=f[k-1][rstar])
{
x[k]=1;
rstar=rstar-w[k];
}
else
{
x[k]=0;
}
}
x[0]=ma;
for(i=0;i<n+1;i++)
{
free(f[i]);
}
free(f);
free(fend);
return x;
}
int *DPQKPMod(int **aintP1, int *aintW, double *aintS0, int intNumofVars, int intKpCapacity,int *aintX)
{
int i=0,j=0,k=0, l=0,r=0;
int intSol=0;
int intBeta=0;
int intExact=0;
double gap=0;
struct timeb start, end;
int *aintFzero = NULL;
int *aintE1 = NULL;
int *aintOpt1=NULL;
int *aintInd=NULL;
int **aintP2Ones=NULL;
int **aintBZeros=NULL;
FILE *MyFile;
// MyFile=fopen("QKP160.txt","a");
aintE1 = (int*)calloc(intNumofVars,sizeof(int));
aintOpt1 = (int*)calloc(intNumofVars,sizeof(int));
aintInd = (int*)calloc(intNumofVars,sizeof(int));
aintP2Ones = (int**)calloc(intNumofVars, sizeof(int*));
aintFzero = (int*)calloc(intKpCapacity+1,sizeof(int));
aintP2Ones = (int**)calloc(intNumofVars, sizeof(int*));
aintBZeros = (int**)calloc(intNumofVars, sizeof(int*));
ftime(&start);
for(i=0; i<intNumofVars; i++)
{
aintP2Ones[i] = (int*)calloc(intNumofVars, sizeof(int));
aintBZeros[i] = (int*)calloc(intKpCapacity+1,sizeof(int));
aintInd[i]=i;
aintX[i] = 0;
for(j=0; j<intNumofVars; j++)
{
aintP2Ones[i][j] = 1;
}
}
aintInd=MySort(aintS0,intNumofVars,aintInd);
int p=0;
for(p=0; p<intNumofVars; p++)
{
k=aintInd[p];
// printf("%d\t",k);
for(r=intKpCapacity; r>-1; r--)
{
if(r >= aintW[k])
{
for(i=0; i<intNumofVars; i++)
{
aintE1[i] = aintBZeros[i][r-aintW[k]];
}
aintE1[k] = 1;
intBeta=aintFzero[r-aintW[k]]+aintP1[k][k];
for(i=0; i<p; i++)
{
if(i!=k && aintBZeros[aintInd[i]][r-aintW[k]]==1)
{
intBeta=intBeta+2*aintP1[k][aintInd[i]];
}
}
int SumonesE1=0;
int SumonesBkr=0;
for(i=0;i<intNumofVars;i++)
{
if(aintE1[i]==1)SumonesE1++;
if(aintBZeros[i][r]==1)SumonesBkr++;
}
//tie-break
if((intBeta > aintFzero[r]) || (intBeta==aintFzero[r] && SumonesE1-1>=SumonesBkr))
{
// printf("ENtry %d\t%d\n",aintFzero[r],intBeta);
aintFzero[r] = intBeta;
for(i=0; i<intNumofVars; i++)
{
aintBZeros[i][r] = aintE1[i];
}
// printf("Exit %d\t%d\n",aintFzero[r],intBeta);
}
}
}
}
intSol = aintFzero[0];
int intrstar=0;
for(i=0; i<intKpCapacity+1; i++)
{
if(intSol<=aintFzero[i])
{
intSol = aintFzero[i];
intrstar=i;
}
}
for(l=0; l<intNumofVars; l++)
{
aintOpt1[l] = aintBZeros[l][intrstar];
}
int Solu=0;
intExact= 0;
for(i=0;i<intNumofVars;i++)
{
for(j=0;j<intNumofVars;j++)
{
intExact += aintX[i]*aintP1[i][j]*aintX[j];
Solu=Solu+((aintP1[i][j]*aintOpt1[i])*aintOpt1[j]);
}
}
ftime(&end);
// printf("Solu=%d\t and Exact=%d the time taken is:%.3f\n",Solu, intExact,(end.time+(double)end.millitm/1000)-(start.time+(double)start.millitm/1000));
int ViolSol=0;
int ViolExac=0;
// fprintf(MyFile,"-----------Optimizer:\n");
for(i=0;i<intNumofVars;i++)
{
ViolSol=ViolSol+aintW[i]*aintOpt1[i];
ViolExac=ViolExac+aintW[i]*aintX[i];
// fprintf(MyFile,"%d\t", aintOpt1[i]);
aintX[i]=aintOpt1[i];
}
/* gap=((double) (intExact-intSol)/intExact)*100;
fprintf(MyFile,"\n Opt Solution: %d \n", intSol);
fprintf(MyFile,"\n Exact Sol: %d \n",intExact);
fprintf(MyFile,"\n gap: %f \n", gap);
printf("Has ours Violate constraints??%d\n",ViolSol);
printf("Has the exact Violate constraints??%d\n",ViolExac);
fprintf(MyFile,"Ellapsed Times: %.3f\n",(end.time+(double)end.millitm/1000)-(start.time+(double)start.millitm/1000));*/
for(i=0; i<intNumofVars; i++)
{
free(aintP2Ones[i]);
free(aintBZeros[i]);
}
free(aintP2Ones);
free(aintBZeros);
free(aintE1);
free(aintOpt1);
free(aintFzero);
free(aintInd);
// fclose(MyFile);
return aintX;
}
int main(void)
{
int Nblock=1;
int intNumofVars=40;
struct timeb starttime, endtime,timestartLB3,timeendLB3, timeFE;
double timeLB3;
int intKpCapacity = 0, SolFE = 0;
int i=0,j=0,k=0,remC3=0;
int LB3 = 0;
double gapLB3= 0.0;
double *S01=NULL;
int *Pi13=NULL;
int *Pi1=NULL;
int *w3=NULL;
int *w1=NULL;
double *B3=NULL;
double *TimeLB3=NULL;
double *GapLB3=NULL;
double *TimeFE=NULL;
double *GapFE=NULL;
int *Ind=NULL;
int *Alt=NULL;
int *Alt3=NULL;
int **aintP1=NULL;
int *aintW=NULL;
double *aintS3=NULL;
double *aintS4=NULL;
int *P3=NULL;
int *P4=NULL;
int *w33=NULL;
int *w44=NULL;
int *X=NULL;
int *X3=NULL;
int *X4=NULL;
int *XFE=NULL;
FILE *input;
FILE *output;
input=fopen("fileread.txt","r");
output=fopen("OutputFile.txt","w");
int b=0;
int Exact=0;
char size[5];
char a[12];
char d[2];
int N;
double t3=0;
B3=(double*)calloc(2,sizeof(double));
TimeLB3=(double*)calloc(Nblock,sizeof(double));
GapLB3=(double*)calloc(Nblock,sizeof(double));
TimeFE=(double*)calloc(Nblock,sizeof(double));
GapFE=(double*)calloc(Nblock,sizeof(double));
S01=(double*)calloc(intNumofVars,sizeof(double));
Pi13=(int*)calloc(intNumofVars,sizeof(int));
Pi1=(int*)calloc(intNumofVars,sizeof(int));
w3=(int*)calloc(intNumofVars,sizeof(int));
w1=(int*)calloc(intNumofVars,sizeof(int));
P3=(int*)calloc(intNumofVars+1,sizeof(int));
w33=(int*)calloc(intNumofVars+1,sizeof(int));
P4=(int*)calloc(intNumofVars+1,sizeof(int));
w44=(int*)calloc(intNumofVars+1,sizeof(int));
X=(int*)calloc(intNumofVars,sizeof(int));
X3=(int*)calloc(intNumofVars,sizeof(int));
X4=(int*)calloc(intNumofVars+1,sizeof(int));
XFE=(int*)calloc(intNumofVars,sizeof(int));
Ind=(int*)calloc(intNumofVars,sizeof(int));
Alt=(int*)calloc(intNumofVars-1,sizeof(int));
Alt3=(int*)calloc(intNumofVars-1,sizeof(int));
aintP1=(int**)calloc(intNumofVars,sizeof(int*));
aintW=(int*)calloc(intNumofVars,sizeof(int));
aintS3=(double*)calloc(intNumofVars,sizeof(double));
aintS4=(double*)calloc(intNumofVars,sizeof(double));
if(input==NULL)
{
printf("error openning the input file");
}
else
{
while(b<Nblock)
{
fscanf(input,"%12s",a);
fscanf(input,"%5s",size);
fscanf(input,"%d",&N);
for(i=0;i<intNumofVars;i++)
{
fscanf(input,"%d",&aintW[i]);
aintP1[i]=(int*)calloc(intNumofVars,sizeof(int));
}
for(i=0;i<intNumofVars;i++)
{
for(j=0;j<intNumofVars;j++)
{
fscanf(input,"%d",&aintP1[i][j]);
}
}
//printf("\n");
for(i=0;i<intNumofVars;i++)
{
fscanf(input,"%d",&X[i]);
// printf(" %d ",X[i]);
}
// printf("\n");
fscanf(input,"%2s",d);
fscanf(input,"%d",&intKpCapacity);
Exact=0;
/********************** \Pi^3_i/w_i *****************************/
#if PI3
ftime(&starttime);
for(i=0;i<intNumofVars;i++)
{
k=0;
remC3=intKpCapacity;
for(j=0;j<intNumofVars;j++)
{
if(j!=i)
{
Alt3[k]=j;
w3[k]=aintW[j];
Pi13[k]=aintP1[i][j];
k++;
}
}
B3=ContinuousKP(Pi13,w3,remC3,intNumofVars-1,B3);
//B3=Dynamic(Pi13,w3,remC3,intNumofVars-1,B3);
aintS3[i]=(aintP1[i][i]+B3[0])/aintW[i];
P3[i+1]=aintP1[i][i]+B3[0];
w33[i+1]=aintW[i];
//printf("%3.f\t%.3f\t%d\n",B3[0],aintS3[i],X[i]);
}
ftime(×tartLB3);
X3=DPQKPMod(aintP1,aintW,aintS3,intNumofVars,intKpCapacity,X3);
ftime(&endtime);
SolFE = 0;
Exact = 0;
LB3=0;
t3=(endtime.time+(double)endtime.millitm/1000)-(starttime.time+(double)starttime.millitm/1000);
fprintf(output,"\n DP solution without fill-up-exchange\n");
for(i=0;i<intNumofVars;i++)
{
fprintf(output," %d ", X3[i]);
for(j=0;j<intNumofVars;j++)
{
LB3 += X3[i]*aintP1[i][j]*X3[j];
}
}
X3=FillUpExchange(X3,aintP1,aintW,intKpCapacity,intNumofVars);
ftime(&timeFE);
fprintf(output,"\n DP solution after fill-up-exchange\n");
for(i=0;i<intNumofVars;i++)
{
fprintf(output," %d ", X3[i]);
for(j=0;j<intNumofVars;j++)
{
SolFE += X3[i]*aintP1[i][j]*X3[j];
Exact+=X[i]*aintP1[i][j]*X[j];
}
}
ftime(&timeendLB3);
timeLB3=(timestartLB3.time+(double)timestartLB3.millitm/1000)-(starttime.time+(double)starttime.millitm/1000)+ (timeendLB3.time+(double)timeendLB3.millitm/1000)-(endtime.time+(double)endtime.millitm/1000);
gapLB3=(double)(Exact-LB3)/Exact*100;
GapFE[b] = (double)(Exact-SolFE)/Exact*100;
TimeFE[b] = (timeFE.time+(double)timeFE.millitm/1000)-(starttime.time+(double)starttime.millitm/1000);
GapLB3[b]=gapLB3;
TimeLB3[b]=timeLB3;
printf("\n\n timeLB3 = %.3f TimeFE =%.3f \n",TimeLB3[b], TimeFE[b]);
printf("\n\n LB3 %d SolFE=%d\t and Exact Now =%d GapFE = %.3f\n",LB3,SolFE,Exact, GapFE[b]);
#endif
/*************** \pi^4 Ordering******************/
#if PI4
ftime(&starttime);
for(i=0;i<intNumofVars;i++)
{
k=0;
remC=intKpCapacity;//-aintW[i];
for(j=0;j<intNumofVars;j++)
{
if(j!=i)
{
Alt[k]=j;
w1[k+1]=aintW[j];
Pi1[k+1]=aintP1[i][j];
k++;
}
}
Ind=Dynamic(Pi1,w1,remC,intNumofVars-1,Ind);
aintS4[i]=(aintP1[i][i]+Ind[0])/aintW[i];
P4[i+1]=aintP1[i][i]+Ind[0];
w44[i+1]=aintW[i];
}
ftime(×tartLB5);
//LB8=DPQKPMod(aintP1,aintW,aintS4,intNumofVars,intKpCapacity,X);
ftime(&endtime);
t2=(endtime.time+(double)endtime.millitm/1000)-(starttime.time+(double)starttime.millitm/1000);
Time8[b]=t2;
printf("Time for pi^2_i=%.5f\n",t2);
X4=Dynamic(P4,w44,intKpCapacity,intNumofVars,X4);
UB4=X4[0];
for(i=0;i<intNumofVars;i++)
{
for(j=0;j<intNumofVars;j++)
{
LB4=LB4+X4[i+1]*aintP1[i][j]*X4[j+1];
}
}
ftime(&timeendLB5);
timeLB5=(timestartLB5.time+(double)timestartLB5.millitm/1000)-(starttime.time+(double)starttime.millitm/1000)+ (timeendLB5.time+(double)timeendLB5.millitm/1000)-(endtime.time+(double)endtime.millitm/1000);
gapUB5=(double)(UB4-Exact)/Exact*100;
gapLB5=(double)(Exact-LB4)/Exact*100;
gap8=(double)(Exact-LB8)/Exact*100;
Gap8[b]=gap8;
GapLB5[b]=gapLB5;
TimeLB5[b]=timeLB5;
printf("Iteration %d LB4=%.2f\t and TimeLB5=%.5f\n",b,LB4,timeLB5);
LB4=0;
#endif
b++;
}
}
free(B3);
free(TimeLB3);
free(GapLB3);
free(TimeFE);
free(GapFE);
free(S01);
free(Pi13);
free(Pi1);
free(w1);
free(w3);
free(P3);
free(w33);
free(P4);
free(w44);
free(X);
free(X3);
free(X4);
free(XFE);
free(Ind);
free(Alt);
free(Alt3);
free(aintW);
free(aintS3);
free(aintS4);
for(i=0;i<intNumofVars;i++)
{
free(aintP1[i]);
}
free(aintP1);
fclose(input);
fclose(output);
return 0;
}