Mahalanobis Indexes for optimal Mahalanobis partition
Loads K partitions and for each of them calculatse indices

(* 
Input: set of K partitions obtained using the MInc[] module (each partition consists of clusters with corresponding weights associated to their elements);
Output: MCH[s], MDB[s], MSSWC[s], and MArea[s] indexes for each
s=2,...,K, with graphs;
Submodules: dM[], MWF[], Hdist[], DIRECT[], MWKMeans[]
*)

Modules

Squared Mahalanobis distance

n=2;
dM[xpod_,ypod_,S_]:=(Det[S]^(1/n))(xpod-ypod).Inverse[S].(xpod-ypod);

Mahalanobis objective function

MWF[PI_,c_,Scov_]:=Module[{m=Length[Flatten[PI,1]]},
Return[Sum[PI[[j,s,1]]*dM[PI[[j,s,2]],c[[j]],Scov[[j]]],{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

Mahalanobis Calinski - Harabasz index (F je Area)

Clear[MCH]
MCH[PI_,z_,Scov_]:=Module[{k=Length[z],F,G},
m=Sum[Length[PI[[j]]],{j,k}];
F=MWF[PI,z,Scov];  (* or *)F=n Sum[Total[PI[[j,All,1]]] (Det[Scov[[j]]])^(1/n),{j,Length[PI]}];
G=Sum[Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}]*dM[z[[j]],cc,Scov[[j]]],{j,k}];
(m-k)*G/((k-1)F)
                                     ]

Mahalanobis Davies - Bouldin index

Clear[MDB]
MDB[PI_,z_,Scov_]:=Module[{k=Length[z],FF,VV,mj,sum=0},
FF=VV=Table[0,{j,k}];
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,i,1]]*dM[PI[[j,i,2]],z[[j]],Scov[[j]]],{i,mj}];
VV[[j]]=FF[[j]]/Sum[PI[[j,i,1]],{i,mj}]
   ,{j,k}];
Do[
max={};
Do[
If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/dM[z[[i]],z[[j]],Scov[[j]]] ]]
,{i,k}];
sum=sum+Max[max]
,{j,k}];
sum/k
            ]

Mahalanobis Silhouette With Criterion

Clear[MSWC]
MSWC[PI_,z_,Scov_]:=Module[{k=Length[PI],alfe,bete,w,s},
w=Table[Total[PI[[i,All,1]]],{i,k}];
alfe=Table[Sum[PI[[i,s,1]]*dM[PI[[i,j,2]],PI[[i,s,2]],Scov[[i]]],{s,Length[PI[[i]]]}]/w[[i]],{i,k},{j,Length[PI[[i]]]}];
bete=Table[Min[Table[Sum[PI[[q,s,1]]*dM[PI[[i,j,2]],PI[[q,s,2]],Scov[[q]]],{s,Length[PI[[q]]]}]/w[[q]],{q,Delete[Range[k],i]}]],{i,k},{j,Length[PI[[i]]]}];
s=Table[(bete[[i,j]]-alfe[[i,j]])/Max[alfe[[i,j]],bete[[i,j]]],{i,k},{j,Length[PI[[i]]]}];
Return[N[Mean[Flatten[s]]]];
]

Mahalanobis Simplify Silhouette With Criterion

MSSWC[PI_,c_,Scov_]:=Module[{k=Length[PI],alfe,bete,s,W},
alfe=Table[dM[PI[[i,j,2]],c[[i]],Scov[[i]]],{i,k},{j,Length[PI[[i]]]}];
bete=Table[Min[Table[dM[PI[[i,j,2]],c[[q]],Scov[[q]]],{q,Delete[Range[k],i]}]],{i,k},{j,Length[PI[[i]]]}];
s=Table[PI[[i,j,1]](bete[[i,j]]-alfe[[i,j]])/Max[alfe[[i,j]],bete[[i,j]]],{i,k},{j,Length[PI[[i]]]}];
(*Return[N[Mean[Flatten[s]]]];*)
W=Total[Table[Total[PI[[i,All,1]]],{i,k}]];
Return[N[Total[Flatten[s]]/W]];
]

Mahalanobis Area index

Clear[MArea]
MArea[PI_,z_,Scov_]:=Module[{},
(*Sum[Det[Scov[[j]]]/Length[PI[[j]]],{j,Length[PI]}]*)
Sum[Sqrt[Det[Scov[[j]]]]/Total[PI[[j,All,1]]],{j,Length[PI]}]
]

Louding of partitions and calculating centers and covariance matrices

SetDirectory[NotebookDirectory[]]; 
PI=ReadList["PART.txt"];

							(* Centers *)
n=Length[PI[[1,1,1,2]]]; 
centri=Table[0,{j,Length[PI]}];
Do[
centri[[br]]=Table[Sum[PI[[br,j,s,1]]*PI[[br,j,s,2]],{s,Length[PI[[br,j]]]}]/Sum[PI[[br,j,s,1]],{s,Length[PI[[br,j]]]}],{j,Length[PI[[br]]]}]
,{br,Length[PI]}]

					(* Covariance Matrices Scov *)
pom=mj=Table[0,{j,Length[PI]+Length[PI[[1]]]-1}];Scov=Table[Table[0,{s,Length[PI[[j]]]}],{j,Length[PI]}];
Do[
Do[
pom[[j]]=Table[Sqrt[PI[[br,j,s,1]]](PI[[br,j,All,2]][[s]]-centri[[br,j]]),{s,Length[PI[[br,j]]]}];
mj[[j]]= Sum[PI[[br,j,s,1]],{s,Length[PI[[br,j]]]}]
,{j,Length[PI[[br]]]}];
Do[Scov[[br,j]]=Table[pom[[j,All,s]].pom[[j,All,r]],{s,n},{r,n}]/mj[[j]]
,{j,Length[PI[[br]]]}];
,{br,Length[PI]}]

							(* Center of the set A *)
Pod=Flatten[PI[[1]],1]; m=Length[Pod];
cc=Total[Pod[[All,1]]*Pod[[All,2]]]/Total[Pod[[All,1]]];

Indexes
MSWC takes a very long time, and the results are similar as by MSSWC

Do[i
Print[i," ",centri[[i]],"; k= ",Length[PI[[i]]],"; mj:",Table[Length[PI[[i,s]]],{s,Length[PI[[i]]]}],"; F=",MWF[PI[[i]],centri[[i]],Scov[[i]]]],
{i,Length[PI]}]

Figures

IndF=IndArea=IndMCH=IndMSSWC=IndSWC=IndMDB={};
br0=Length[PI[[1]]]-1;
time=Timing[Do[
IndF=Append[IndF,MWF[PI[[br]],centri[[br]],Scov[[br]]]]; 
,{br,Length[PI]}]]; length=Length[IndF];
Print["IndF: ",TabIndF=Table[{j+br0,IndF[[j]]},{j,length}],";  Time=",time[[1]]," sec"];
time=Timing[Do[
IndArea=Append[IndArea,MArea[PI[[br]],centri[[br]],Scov[[br]]]];   
,{br,Length[PI]}]];
Print["IndArea: ",TabIndArea=Table[{j+br0,IndArea[[j]]},{j,length}],";  Time=",time[[1]]," sec"];
time=Timing[Do[
IndMCH=Append[IndMCH,MCH[PI[[br]],centri[[br]],Scov[[br]]]];   
,{br,Length[PI]}]];
Print["IndMCH: ",TabIndMCH=Table[{j+br0,IndMCH[[j]]},{j,length}],";  Time=",time[[1]]," sec"];
time=Timing[Do[
IndMDB=Append[IndMDB,MDB[PI[[br]],centri[[br]],Scov[[br]]]];   
,{br,Length[PI]}]];
Print["IndMDB: ",TabIndMDB=Table[{j+br0,IndMDB[[j]]},{j,length}],";  Time=",time[[1]]," sec"];
time=Timing[Do[
IndMSSWC=Append[IndMSSWC,MSSWC[PI[[br]],centri[[br]],Scov[[br]]]];    
,{br,Length[PI]}]];
Print["IndMSSWC: ",TabIndMSSWC=Table[{j+br0,IndMSSWC[[j]]},{j,length}],";  Time=",time[[1]]," sec"];
(*
time=Timing[Do[
IndSWC=Append[IndSWC,SWC[PI[[br]],centri[[br]],Scov[[br]]]];     
,{br,Length[PI]}]];
Print["IndSWC: ",TabIndSWC=Table[{j+br0,IndSWC[[j]]},{j,length}],";  Time=",time[[1]]," sec"];
*)
									(* Slike *)
slF=ListPlot[TabIndF,Joined->True,AxesOrigin->{1,.75Min[IndF]},PlotLabel->"F",ImageSize->200];
slArea=ListPlot[TabIndArea,Joined->True,AxesOrigin->{1,.75Min[IndArea]},PlotLabel->"MArea",ImageSize->200];
slMCH=ListPlot[TabIndMCH,Joined->True,AxesOrigin->{1,.75Min[IndMCH]},PlotLabel->"MCH",ImageSize->200];
slMDB=ListPlot[TabIndMDB,Joined->True,AxesOrigin->{1,.75Min[IndMDB]},PlotLabel->"MDB",ImageSize->200];
slMSSWC=ListPlot[TabIndMSSWC,Joined->True,AxesOrigin->{1,.75Min[IndMSSWC]},PlotLabel->"MSSWC",ImageSize->200];
Print[Style[GraphicsGrid[{{slF,slMCH,slMDB},{slMSSWC,slArea}}],Magnification->1.2]]

Figures

slArea=ListPlot[TabIndArea,Joined->True,AxesOrigin->{1,.75Min[IndArea]},PlotLabel->"MArea",ImageSize->200];
slMCH=ListPlot[TabIndMCH,Joined->True,AxesOrigin->{1,.75Min[IndMCH]},PlotLabel->"MCH",ImageSize->200];
slMDB=ListPlot[TabIndMDB,Joined->True,AxesOrigin->{1,.75Min[IndMDB]},PlotLabel->"MDB",ImageSize->200];
slMSSWC=ListPlot[TabIndMSSWC,Joined->True,AxesOrigin->{1,.75Min[IndMSSWC]},PlotLabel->"MSSWC",ImageSize->200];
Print[Style[GraphicsGrid[{{slMCH,slMDB,slMSSWC,slArea}}],Magnification->1.2]]

Figures corrected

slArea=ListPlot[TabIndArea,Joined->True,AxesOrigin->{1,0},Ticks->{Automatic,{.01,.02,.03}},GridLines->{Automatic,{0,.01,.02,.03,.04}},PlotRange->{0,.03}];
slVCH=ListPlot[TabIndMCH,Joined->True,AxesOrigin->{1,1000},Ticks->{Automatic,{1000,2000,3000,4000,5000}},GridLines->{Automatic,{1000,2000,3000,4000,5000}},PlotRange->{1000,5000}];
slVDB=ListPlot[TabIndMDB,Joined->True,AxesOrigin->{1,0},Ticks->{Automatic,{.1,.2,.3,.4}},GridLines->{Automatic,{.1,.2,.3,.4}},PlotRange->{0,.4}];
slSWC=ListPlot[TabIndMSSWC,Joined->True,AxesOrigin->{1,.8},GridLines->{Automatic,{.825,.85,.875,.90}},PlotRange->{.8,.9},Ticks->{Automatic,{.75,.8,.85,.90}}];
Print[GraphicsGrid[{{slVCH,slVDB,slSWC,slArea}},ImageSize->Large]]