
#############################################################
#  procedure name: BWM-v1 by Huseyin Kocak 03/03/2017       #
#  last modified-v1.1: total deviation added 17/07/2017     #
#############################################################
#  BO represents best-to-others comparison vector           #
#  OW represents others-to-worst comparison vector          #
#############################################################


BWM:=proc(BO,OW)

local n,n2,abest,abest2,aworst,aworst2,nn1,nn2,nn3,nn4,i,best1,worst1,best2,worst2,sumBO,sumOW,sumBO2,sumOW2,sumw,sumw1,sol,CI,bob,bow,owb,oww,nn5,nn6,nn7,nn8,sumBOL,sumOWL,sumBO2L,sumOW2L,sol1,CI1,TD1,TD2,sumTD11,sumTD12,sumTD21,sumTD22;


#BO:=readstat("please enter your BO vector as [.,.,.,.]");
#OW:=readstat("please enter your OW vector as [.,.,.,.]");

## check the sizes of the vectors
n:=size(BO,2);
n2:=size(OW,2);
if abs(n-n2)>0 then error "sizes of BO and OW should be same and equal to number of criteria"; end if;
print("the number of the criteria",n);

## check if a_BB=1, a_WW=1 and value of a_BW same in both vector
abest:=min(BO);
aworst:=max(BO);
abest2:=max(OW);
aworst2:=min(OW);
if abest<>1 then error "best-best comparison should be 1"; end if;
if aworst2<>1 then error "worst-worst comparison should be 1"; end if;
if aworst<>abest2 then error "value of the best-worst comparison should be same in BO and OW"; end if;

## be sure that the number of best and worst criteria is same in both vector
## also be sure that best and worst criteria are same in both vector
nn1:=0;nn2:=0;nn3:=0;nn4:=0;
for i from 1 to n do
if BO[i]=abest then nn1:=nn1+1; bob[nn1]:=i; end if;
if BO[i]=aworst then nn2:=nn2+1; bow[nn2]:=i; end if;
if OW[i]=abest2 then nn3:=nn3+1; owb[nn3]:=i; end if;
if OW[i]=aworst2 then nn4:=nn4+1; oww[nn4]:=i; end if;
end do;
if nn1<>nn3 then error "number of best should be same in BO and OW"; end if;
if nn2<>nn4 then error "number of worst should be same in BO and OW"; end if;
if bob[nn1]<>owb[nn3] then error "best criteria are not same in BO and OW"; end if;
if bow[nn2]<>oww[nn4] then error "worst criteria are not same in BO and OW"; end if;


if nn1=1 then print("the best criterion",C[bob[nn1]]); else 
print("the best criteria",seq(C[bob[j]],j=1..nn1));
end if;
if nn2=1 then print("the worst criterion",C[bow[nn2]]); else 
print("the worst criteria",seq(C[bow[j]],j=1..nn2));
end if;
print("the value of the best-worst comparison",aworst);

## Euclidean BWM computation
sumBO:=0;
sumOW:=0;
sumBO2:=0;
sumOW2:=0;
sumw:=0;
sumw1:=0;
for i from 1 to n do

if i<>bob[nn1] then sumBO:=sumBO+(w1[bob[nn1]]/w1[i]-BO[i])^2;
end if;
if i<>bob[nn1] then if i<>bow[nn2] then sumOW:=sumOW+(w1[i]/w1[bow[nn2]]-OW[i])^2; end if; end if;
sumw:=sumw+w1[i];

if i<>bob[nn1] then sumBO2:=sumBO2+(w2[bob[nn1]]/w2[i]-aworst)^2;
end if;
if i<>bob[nn1] then if i<>bow[nn2] then sumOW2:=sumOW2+(w2[i]/w2[bow[nn2]]-aworst)^2; end if; end if;
sumw1:=sumw1+w2[i];

end do;

sol:=Minimize(sumBO+sumOW,{sumw=1}, assume = nonnegative);
assign(sol[2]);
CI:=Minimize(sumBO2+sumOW2,{sumw1=1}, assume = nonnegative);

#TD calculation for Euclidean
sumTD11:=0;
sumTD12:=0;
for i from 1 to n do
sumTD11:=sumTD11+(BO[i]-(w1[bob[nn1]]/w1[i]))^2;
sumTD12:=sumTD12+(OW[i]-(w1[i]/w1[bow[nn2]]))^2;
end do;

print("the Euclidean BWM results are as follows:");
for i from 1 to n do
printf("\t\r w%a=%a",i,w1[i]);
end do;
printf("\t\r xi=%a",sqrt(sol[1]));
printf("\t\r consistency index=%a",sqrt(CI[1]));
printf("\t\r the consistency ratio=xi/consistency index=%a",sqrt(sol[1]/CI[1]));
printf("\t\r Total Deviation=%a",(sumTD11+sumTD12)/(2*n));

## linear Chebyshev BWM computation

sumw:=0;
sumw1:=0;
nn5:=0;nn6:=0;nn7:=0;nn8:=0;
for i from 1 to n do

if i<>bob[nn1] then
nn5:=nn5+2;
sumBOL[nn5-1]:=(w3[bob[nn1]]-BO[i]*w3[i])<=x;
sumBOL[nn5]:=(w3[bob[nn1]]-BO[i]*w3[i])>=-x;
end if;
if i<>bob[nn1] then if i<>bow[nn2] then
nn6:=nn6+2;
sumOWL[nn6-1]:=(w3[i]-OW[i]*w3[bow[nn2]])<=x;
sumOWL[nn6]:=(w3[i]-OW[i]*w3[bow[nn2]])>=-x;
end if; end if;
sumw:=sumw+w3[i];

if i<>bob[nn1] then 
nn7:=nn7+2;
sumBO2L[nn7-1]:=(w4[bob[nn1]]-aworst*w4[i])<=y;
sumBO2L[nn7]:=(w4[bob[nn1]]-aworst*w4[i])>=-y;
end if;
if i<>bob[nn1] then if i<>bow[nn2] then
nn8:=nn8+2;
sumOW2L[nn8-1]:=(w4[i]-aworst*w4[bow[nn2]])<=y;
sumOW2L[nn8]:=(w4[i]-aworst*w4[bow[nn2]])>=-y;
end if; end if;
sumw1:=sumw1+w4[i];

end do;

sol1:=LPSolve(x,{seq(sumBOL[j],j=1..nn5),seq(sumOWL[j],j=1..nn6),sumw=1,seq(w3[j]>=0,j=1..n)});
assign(sol1[2]);
CI1:=LPSolve(y,{seq(sumBO2L[j],j=1..nn7),seq(sumOW2L[j],j=1..nn8),sumw1=1,seq(w4[j]>=0,j=1..n)});

#TD calculation for Chebyshev
sumTD21:=0;
sumTD22:=0;
for i from 1 to n do
sumTD21:=sumTD21+(BO[i]-(w3[bob[nn1]]/w3[i]))^2;
sumTD22:=sumTD22+(OW[i]-(w3[i]/w3[bow[nn2]]))^2;
end do;

print("the linear Chebyshev BWM results are as follows:");
for i from 1 to n do
printf("\t\r w%a=%a",i,w3[i]);
end do;
printf("\t\r xi=%a",sol1[1]);
printf("\t\r max xi=%a",(CI1[1]));
printf("\t\r the actual consistency ratio=xi/max xi=%a",(sol1[1]/CI1[1]));
printf("\t\r Total Deviation=%a",(sumTD21+sumTD22)/(2*n));

end proc:
