The program for the manuscript (FuzzyInc-New): 
  A partitioning incremental algorithm  using adaptive Mahalanobis fuzzy clustering

Fuzzy Mahalanobis Clustering

(* Fuzzy Mahalanobis Clustering is a complete Mathematica-program for
fuzzy Mahalanobis clustering using the normalized Mahalanobis distance-like function.All modules needed are listed in the preamble.One has
to load one of the suggested examples first.Then one can perform the
GKc-means algorithm GKcmeans[] or the Mahalanobis fuzzy incremental
algorithm MFInc[],which will,based on several fuzzy indexes,propose one
of the partitions obtained as the MAPart.If the original partition is known,it will also calculate the Rand and Jaccard indexes,the confusion matrix,and the Hausdorff distance between sets of original and calculated cluster
centers.
*)

Modules

M - distance and Membership matrix

dm[x_,y_,S_]:=(x-y).Inverse[S].(x-y)
dM[x_,y_,S_]:=(Det[S]^(1/n))(x-y).Inverse[S].(x-y)

muGK[A_,cen_,Scov_,i_,j_]:=1/(Sum[(dM[A[[i]],cen[[j]],Scov[[j]]]/dM[A[[i]],cen[[s]],Scov[[s]]])^(1/(q-1)),{s,Length[cen]}])

Fuzzy Mahalanobis objective function

MFz[A_,cen_,U_,cov_]:=Module[{m=Length[A],k=Length[cen]},
Return[
Sum[muGK[A,cen,cov,i,j]^q dM[cen[[j]],A[[i]],cov[[j]]],{i,m},{j,k}]
	]
];
MF[A_,cen_,U_,cov_]:=Module[{m=Length[A],k=Length[cen]},
Return[
Sum[Min[Table[muGK[A,cen,cov,i,j]^q dM[cen[[j]],A[[i]],cov[[j]]],{j,k}]]
,{i,m}]
	]
];

Mahalanobis objective function

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

Minimal element of the upper triangular matrix

Input: The upper triangular matrix ℝ;
Output: Minimal element;
Submodules: None

MinR[R_]:=Module[{n=Length[R[[1]]],min=Infinity},
Do[
red=Drop[R[[i]],{1,i}];
minn=Min[red];
If[minn<min,min=minn];
,{i,n-1}];
min]

Testing if two ellipsoids is not intersecting

Input: Ell1[c1,cov1],Ell2[c2,cov2];
Output: True (not intersecting)or False (intersecting)
Submodules: None

IntersectTest[c1_,cov1_,c2_,cov2_]:=Module[{FF,max},
FF[lambda_]:=lambda (1-lambda) (c2-c1).Inverse[lambda cov2+(1-lambda) cov1].(c2-c1);
max=First[FindMaximum[{FF[x],0<x<1},x]];
max>1
]

GKcmeans modul

GKcmeans[A_,z_,Ind_]:=Module[{m=Length[A],k=Length[z],U2,slcen,n=Length[A[[1]]],eps=0.005,It=20},  
If[Ind≠0,slA=ListPlot[A]];
					(* Initial approximation *)
kk=1;
cen=z//N; 
Scov=ScovN=Table[IdentityMatrix[n],{j,k}]; 
U1=Table[muGK[A,cen,Scov,i,j],{i,m},{j,k}];  
If[n==2,slcen=ListPlot[cen,PlotStyle->{Red,PointSize[.04],Opacity[.7]}];
sl2=Show[slA,slcen,PlotRange->AG,AspectRatio->Automatic,Frame->True,GridLines->Automatic];
Print[GraphicsGrid[{{sl0,sl2}},ImageSize->Medium]]];
F1=  MFz[A,cen,U1,Scov];Print["F0 = ",F1];
						(* While loop *)
While[
cenN=Table[Sum[A[[i]]muGK[A,cen,Scov,i,j]^q ,{i,m}]/Sum[muGK[A,cen,Scov,i,j]^q,{i,m}],{j,k}];
If[Ind≠0,slcen=ListPlot[cenN,PlotStyle->{Red,PointSize[.04],Opacity[.7]}];
sl1=sl2; sl2=Show[slA,slcen,PlotRange->AG,AspectRatio->Automatic,Frame->True,GridLines->Automatic];
Print[GraphicsGrid[{{sl1,sl2}},ImageSize->Medium]]];
Do[ScovN[[j]]=Sum[muGK[A,cenN,Scov,i,j]^q KroneckerProduct[(A[[i]]-cenN[[j]]),(A[[i]]-cenN[[j]])],{i,m}]/Sum[muGK[A,cenN,Scov,i,j]^q,{i,m}]
,{j,k}];
U2=Table[muGK[A,cenN,ScovN,i,j],{i,m},{j,k}]; DelU=U2-U1;
F2=  MFz[A,cenN,U2,ScovN];delF=(F1-F2)/F2;
U1=U2;F1=F2;
If[Ind≠0,Print["It-",kk,": centri=",cenN,"\n||Del U||=",Norm[DelU],"; delF = ",delF ]];
Abs[delF]>eps&& kk<It,
kk=kk+1; cen=cenN; Scov=ScovN;
 ];
{cen,U2,Scov,kk}]

Mahalanobis Fuzzy Incremental algorithm

MFInc[A_,z_,KK_,Ind_:0]:=Module[{k,centri,m=Length[A],n=Length[A[[1]]],dom,TTdir=0,TTGK=0},
					(* Initialization *)
Scov=Table[IdentityMatrix[n],{i,Length[z]}]; 
brk0=Length[cen];
sl2=slA;
FF={};  
centri=z;
ScovN={IdentityMatrix[n]};
PART={}; IndFXB=IndFCH=IndFDB=IndFHV=IndRand=IndJaccard=IndH=IndMFNew={};
				(* Iterative Procedure *)
Do[
Print["\nNa postojecih ",Length[centri]," centara dodaje se 1 novi"];
Clear[u,v,F2]; 
(* Only a few iteration of DIRECT is sufficient because this is a initial approximation for WKMeans *)
					(* Compile *)
FCompile=Compile[{{u,_Real},{v,_Real}},
Evaluate[Sum[Min[Append[Table[muGK[A,centri,ScovN,i,s]dM[A[[i]],centri[[s]],ScovN[[s]]],{s,Length[centri]}],Norm[A[[i]]-{u,v}]^2]],{i,m}] ]];	
F2[u_,v_]:=Apply[FCompile,{u,v}];
								(* DIRECT *)
dom=AG;
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];
cen=Append[centri,min[[1]]];
If[Ind≠0,Print[Length[cen]-1," old centers and one new optimized center: ",Table[cen[[j]],{j,Length[cen]}]]];
k=Length[cen]; 
			(* Call of Adaptive Mahalanobis k-means *)
t1=Timing[solGK=GKcmeans[A,cen,0]];
TTdir=TTdir+minp[[1]];TTGK=TTGK+t1[[1]];
                                            (*  *)
cenN=solGK[[1]];
UN=solGK[[2]];
ScovN=solGK[[3]];
FN=MFz[A,cenN,UN,ScovN];
Print["New Fuzzy GKMeans refined centers: ",cenN];
Print["GKMeans-it = ",solGK[[4]],"\nMFz = ",FN,"; MF = ",MF[A,cenN,UN,ScovN],
"; Time (DIRECT+GKMeans) = ",minp[[1]],"+",t1[[1]]," = ",minp[[1]]+t1[[1]]];
				(* Slika klastera *)
slA=ListPlot[A,PlotStyle->{Brown,Opacity[.4]}];
slPI0=Show[slA,Table[Graphics[{Hue[j/k],Opacity[0.4],Ellipsoid[solGK[[1,j]], kappa solGK[[3,j]]]}],{j,Length[solGK[[1]]]}],GridLines->Automatic,AspectRatio->Automatic,PlotLabel->"Polazna pozicija" ];
					(* Indexes *)
fxb=MFXB[A,solGK[[1]],solGK[[2]],solGK[[3]]];IndFXB=Append[IndFXB,fxb];
fch=MFCH[A,solGK[[1]],solGK[[2]],solGK[[3]]];IndFCH=Append[IndFCH,fch];
fdb=MFDB[A,solGK[[1]],solGK[[2]],solGK[[3]]];IndFDB=Append[IndFDB,fdb];
fhv=MFHV[solGK[[3]]];IndFHV=Append[IndFHV,fhv];
fnew=MFNew1[A,solGK[[1]],solGK[[2]],solGK[[3]],1];IndMFNew=Append[IndMFNew,fnew];
Print["k=",Length[cenN],"; FXB=",fxb,"; FCH=",fch,"; FDB=",fdb,"; FHV=",fhv,"; MFNew=",fnew];
ind=RandFrigue[UN,U];
HD=Hdist[centers,cenN];
IndRand=Append[IndRand,ind[[1]]];
IndJaccard=Append[IndJaccard,ind[[2]]];
IndH=Append[IndH,HD];
If[Ind≠0,Print["Rand: ",IndRand,"\nJaccard: ",IndJaccard,"\nHausdorff: ",IndH]];
				  (* Figures *)  
If[n==2,
k=Length[solGK[[1]]];
slcenN=ListPlot[cenN,PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.04],Opacity[.6]}];
Clear[x,y];  (* This M-circle contains about 68% data points *)
				(* Ovo je samo za crtanje sljedeceg centra *)
					(* Compile *)
Clear[u,v,F2]; 
FCompile=Compile[{{u,_Real},{v,_Real}},
Evaluate[Sum[Min[Append[Table[Norm[A[[i]]-cenN[[s]]]^2,{s,Length[cenN]}],Norm[A[[i]]-{u,v}]^2]],{i,m}] ]];	
F2[u_,v_]:=Apply[FCompile,{u,v}];
mini=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)];
nextc=ListPlot[{mini[[1]]},PlotStyle->{Red,PointSize[.03],Opacity[.8]}];
								(* Figure *)
cl=Table[solGK[[2,All,j]],{j,k}];
boje=Table[(2j-1)/10,{j,k}];
						(* Nijanse *)
 nijanse=Table[3,{j,k}];
slrez=Show[Table[Table[ListPlot[{A[[i]]},PlotStyle->{Hue[boje[[j]]],Opacity[cl[[j,i]]/nijanse[[j]]]}],{i,m}],{j,k}],slcenN,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[slrez];
  ];				 
centri=solGK[[1]];     
PART=Append[PART,{cenN,UN,ScovN,FN}]; 
,{brk,KK}];
								(* Graphs of Indexes *)
Print[Style["prag = ",Blue,14],prag,Style[";  p = ",Blue,14],p];
Print["IndFXB: ",IndFXB,"\nIndFCH: ",IndFCH,"\nIndFDB: ",IndFDB,"\nIndFHV: ",IndFHV,"\nIndMFNew: ",IndMFNew];
Print[GraphicsGrid[{{ListPlot[Table[{i+1,IndFXB[[i]]},{i,Length[IndFXB]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],ListPlot[Table[{i+1,IndFCH[[i]]},{i,Length[IndFCH]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],ListPlot[Table[{i+1,IndFDB[[i]]},{i,Length[IndFDB]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],ListPlot[Table[{i+1,IndFHV[[i]]},{i,Length[IndFHV]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],ListPlot[Table[{i+1,IndMFNew[[i]]},{i,Length[IndMFNew]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic]}},ImageSize->Large]];
Print["RandFrigue: ",IndRand,"\nJaccardFrigue: ",IndJaccard,"\nHausdorff: ",IndH];
{PART,{IndFXB,IndFCH,IndFDB,IndFHV,IndMFNew,IndRand,IndJaccard,IndH},{TTdir,TTGK}}
]

Mahalanobis Fuzzy Indexes

(* Xie-Beni fuzzy index *)
Clear[MFXB]
MFXB[A_,cen_,U_,cov_]:=Module[{m,k},m=Length[A];k=Length[cen];Sum[U[[i,j]]^q dM[cen[[j]],A[[i]],cov[[j]]],{i,m},{j,k}]/(m*Min[Flatten[Table[dM[cen[[i]],cen[[j]],cov[[i]]+cov[[j]]],{i,k-1},{j,i+1,k}],1]])
]

(* Calinski-Harabasz fuzzy index *)
Clear[MFCH]
MFCH[A_,cen_,U_,cov_]:=Module[{m,k},m=Length[A];k=Length[cen];
cc=Mean[A];
covcc=Sum[KroneckerProduct[(A[[i]]-cc),(A[[i]]-cc)],{i,m}]/m;
Return[(Sum[Sum[U[[ii,j]]^q,{ii,m}] *dM[cc,cen[[j]],(*cov[[j]]+*)covcc],{j,k}]/(k-1))/(Sum[U[[i,j]]^q dM[cen[[j]],A[[i]],cov[[j]]],{i,m},{j,k}]/(m-k))]
]

(* Davies-Bouldin fuzzy index *)
Clear[MFDB]
MFDB[A_,cen_,U_,cov_]:=Module[{m=Length[A],k=Length[cen],VV,max,sum=0},
VV=Table[0,{j,k}];
Do[
VV[[j]]=Sum[U[[i,j]]^q dM[cen[[j]],A[[i]],cov[[j]]],{i,m}]/Sum[U[[i,j]]^q,{i,m}]
   ,{j,k}];
Do[
max={};
Do[
If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/dM[cen[[i]],cen[[j]],cov[[i]]+cov[[j]]] ]]
,{i,k}];
sum=sum+Max[max]
,{j,k}];
sum/k
            ]

(* Fuzzy hypervolume index *)
Clear[MFHV]
MFHV[cov_]:=Module[{k=Length[cov],VV,max,sum=0},
Return[Sum[Sqrt[Det[cov[[j]]]],{j,k}]]
]

Fuzzy MInimal Distence index for ellipsoidal clusters

Input: The set 𝒜 ⊂ (ℝ^n), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), Membership matrix U,
covariance matrices (∑_j)∈(ℝ^(n\[Cross]n)), j=1⃛k;
Output: FMMD -index;

Submodules:  MMD[],IntersectTest[],MinR[]
The variable prag∈[0,1]  comes from MFINC[].

(* If u_{ij} ≥ prag, element (Subscript[a, i]) stays in cluster (Subscript[π, j]) *)
MFNew1[A_,Cen_,U_,Cov_,Ind_:0]:=Module[{m=Length[A],k=Length[Cen],n=Length[A[[1]]]},
PI=Table[{},{j,k}]; Res={};
Do[
max=Ordering[U[[i]],-1]; 
If[U[[i,max[[1]]]]>prag,
PI[[max[[1]]]]=Append[PI[[max[[1]]]],A[[i]]],
Res=Append[Res,A[[i]]]
]
,{i,m}];
						(* New centers *)
Cen2=Table[Mean[PI[[j]]],{j,k}];
Print["|PI[[j]]|: ",Table[Length[PI[[j]]],{j,k}],", |Res|=",Length[Res]];
				(* New covariance matrices *)
Cov2=Table[0,{j,k}];
Do[
Cov2[[j]]= Sum[KroneckerProduct[(PI[[j,i]]-Cen2[[j]]),(PI[[j,i]]-Cen2[[j]])],{i,Length[PI[[j]]]}]/Length[PI[[j]]]
,{j,k}];
					(* Call MMD index *)
MMD[PI,Cen2,Cov2,Ind]
]

MInimal Distence index for ellipsoidal clusters

Input: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} of the set 𝒜 ⊂ (ℝ^n), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), 
covariance matrices (∑_j)∈(ℝ^(n×n)), j=1⃛k;
Output: Distance matrix between clusters and minimal distance;
Submodules:  IntersectTest[],MinR[]

Clear[MMD]
MMD[PI_,z_,Scov_,Ind2_:0]:=Module[{k=Length[PI],min,P1,P2,m},
DDD=Table[0,{r,k},{s,k}];
Ell=Table[Ellipsoid[z[[j]],kappa  Scov[[j]]],{j,k}];
slPI=Show[slA,Table[Graphics[{Hue[j/k],Opacity[0.4],Ell[[j]]}],{j,Length[z]}],GridLines->Automatic,AspectRatio->Automatic,PlotLabel->"Novi Cen, Cov, kappa"];
Print[GraphicsGrid[{{slPI0,slPI}},ImageSize->Medium]];
	 (* M-udaljenosti (najprije se testira je li presjek prazan) *)
Do[
Do[
If[IntersectTest[Ell[[r,1]],Ell[[r,2]],Ell[[s,1]],Ell[[s,2]]],
min=FindMinimum[{EuclideanDistance[x,y],{x∈Ell[[r]],y∈Ell[[s]]}}, {x,y}];{P1,P2}={x/.min[[2]],y/.min[[2]]};
Print[{s,r},{{P1,P2},min}];
DDD[[r,s]]=min[[1]],
{P1,P2}={Ell[[r,1]],Ell[[r,1]]};DDD[[r,s]]=0
	];
If[Ind2≠ 0, 
If[n==2,
Print[Show[
Graphics[{boja[[r]],EdgeForm[Directive[Thickness[.004],boja[[r]]]],Opacity[0.075],Ell[[r]]}],
Graphics[{boja[[s]],EdgeForm[Directive[Thickness[.004],boja[[s]]]],Opacity[0.08],Ell[[s]]}],
ListPlot[PI,PlotStyle->Table[{PointSize[.01],boja[[j]]},{j,Length[z]}]],
ListPlot[{P1,P2}, PlotStyle->{Red,PointSize[.025]}], 
Graphics[{Red,Thin,Line[{P1,P2}]}],
Axes->True, AxesOrigin->{0,0},AxesStyle->{Thickness[.005],Thickness[.005]},GridLines->Automatic,PlotRange->AG,AspectRatio->Automatic,ImageSize->Small]],
Print[Show[
Graphics3D[{RGBColor[z[[r]]],Opacity[0.5],Ell[[r]]}],Graphics3D[{RGBColor[z[[s]]],Opacity[0.5],Ell[[s]]}],
ListPointPlot3D[PI0,PlotStyle->Table[{PointSize[.001],Opacity[.5],RGBColor[z[[j]]]},{j,Length[z]}]],
ListPointPlot3D[{P1,P2}, PlotStyle->{Red,PointSize[Large]}],
Graphics3D[{Red,Thin,Line[{P1,P2}]}],
Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},ImageSize->Small]]
		];
			];
,{s,r+1,k}]
,{r,k}];
Print["DDD=",MatrixForm[DDD]," - distance matrix between clusters"];
m=Sum[Length[PI[[j]]],{j,Length[PI]}];
mwf=MWF[PI,z,Scov];
(m-k)MinR[DDD] /mwf
]

Hausdorff distance between two sets of centroids

(* Po defaultu koristi se euklidska metrika. Ako se želi druga, treba to na dva mjesta upisati *)

Hdist[A_,B_]:=Module[{m=Length[A],n=Length[B],u,v},
u=Max[Table[Min[Table[Norm[x-y],{y,B}] ],{x,A}]];
v=Max[Table[Min[Table[Norm[x-y],{x,A}] ],{y,B}]];
Max[u,v]
]

Fuzzy Rand and Jaccard indexes

(* 
Input: membership matrices of partitions (Π^((1))) and (Π^((2))) of the set 𝒜 ⊂ (ℝ^n); 
Output: fuzzy Rand FR and fuzzy Jaccard FJ indexes;
Submodules: none needed 
*)

RandFrigue[A_,B_]:=Module[{m=Dimensions[A][[1]],PsiA,PsiB,a,b,c,d},
PsiA=PsiB=Table[0,{i,m},{j,m}];
Do[
PsiA[[i,j]]=A[[i]].A[[j]];PsiB[[i,j]]=B[[i]].B[[j]];
,{i,m-1},{j,i+1,m}];
a=Sum[PsiA[[i,j]]*PsiB[[i,j]],{i,m-1},{j,i+1,m}];
b=Sum[PsiA[[i,j]]*(1-PsiB[[i,j]]),{i,m-1},{j,i+1,m}];
c=Sum[(1-PsiA[[i,j]])*PsiB[[i,j]],{i,m-1},{j,i+1,m}];
d=Sum[(1-PsiA[[i,j]])*(1-PsiB[[i,j]]),{i,m-1},{j,i+1,m}];
{(a+d)/(a+b+c+d),a/(a+b+c)}
]

Fuzzy index according to Hullermeier et al. (2012)

(* A, B su  membership matrices *)
Hullermeier[A_,B_]:=Module[{m=Dimensions[A][[1]],PsiA,PsiB,dAB},
PsiA=PsiB=Table[0,{i,m},{j,m}];
Do[
PsiA[[i,j]]=1-Norm[A[[i]]-A[[j]],Infinity];PsiB[[i,j]]=1-Norm[B[[i]]-B[[j]],Infinity];
,{i,m-1},{j,i+1,m}];
dAB=Sum[Abs[PsiA[[i,j]]-PsiB[[i,j]]],{i,m-1},{j,i+1,m}]/(m(m-1)/2);
1-dAB
]

Similarity of two ellipses (can be applied on ellipsoids too)

Scn[E1_,E2_]:=Module[{x,y,eig1,eig2,sp,theta,al,be,aa,osi},
eig1=Eigensystem[E1[[3]]];  
eig2=Eigensystem[E2[[3]]]; 
sp=Dot[eig1[[2,1]],eig2[[2,1]]];
theta=ArcCos[sp];
al=Sqrt[eig1[[1]]];be=Sqrt[eig2[[1]]];
(*al=eig1[[1,1]]/eig1[[1,2]];be=eig2[[1,1]]/eig2[[1,2]];*)
osi=Max[Abs[al-be]/al];
aa={Exp[- Sqrt[dM[E1[[1]],E2[[1]],E1[[3]]+E2[[3]]]]], Exp[-Sqrt[2]Abs[Sin[theta]]], Exp[-osi]};
bb=Product[aa[[i]],{i,3}]^(1/3);
{aa,bb}
]

DIRECT Modul

Options[DIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
DIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,ind,novic,novid,novifc,red,TMin,iter=0,dtest,anim={}},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
DIRECT::maxiter="Dosegnut maksimalan broj iteracija `1`.";
n=Length[dom];
(* skaliranje funkcije na hiperkocku duljine stranice 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* početni centar i promjeri *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};
TMin={c[[1]],fc[[1]]};
(* GLAVNA PETLJA *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* KORAK 1: traženje potencijalno optimalnih hiperkvadara *)
dmaxrazni=Union[dmax];
(* svi podaci grupirani po d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* podaci s najmanjim f grupirani po d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* provjera ostalih uvjeta *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* KORAK 2: dijeljenje potencijalno optimalnih *)
Do[
(* smjerovi po kojima se dijeli *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
If[novifc[[2red[[1]]-1]]<TMin[[2]],TMin=N[{novic[[2red[[1]]-1]],novifc[[2red[[1]]-1]]},10]];
If[novifc[[2red[[1]]]]<TMin[[2]],TMin=N[{novic[[2red[[1]]]],novifc[[2red[[1]]]]},10]];
,{pot,pot2}];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* slika za 2D slučaj *)
(*If[n==2,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}]
(* vraćanje u početnu domenu *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* ispis slika za 2D slučaj *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

SymDIRECT modul

(*Za zadani centar i duljine stranica hiperpravokutnika funkcija test vraća:
1 ako se hiperpravokutnik cijel nalazi u području simetrije
0 ako se hiperpravokutnik djelomično nalazi u području simetrije
-1 ako se hiperpravokutnik ne nalazi u području simetrije (ili ako je dodirna samo jedna točka)*)
test[c_,dd_]:=Module[{n=Length[c],S,d=dd/2,unutra=1},
S=Table[1,{n}];
Do[
If[c[[i]]-d[[i]]<c[[i+1]]+d[[i+1]],
unutra=0;Break[];
]
,{i,1,n-1}];
If[unutra==0,
Do[
If[c[[i]]+S[[i]]*d[[i]]<c[[i+1]]+d[[i+1]],
If[c[[i]]+S[[i]]*d[[i]]>c[[i+1]]-d[[i+1]],
S[[i+1]]=-1,
unutra=-1;Break[];
]
]
,{i,1,n-1}];
];
Return[unutra];
]

Options[SymDIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
SymDIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,temp2,ind,novic,novid,novifc,noviunutra,red,TMin,iter=0,dtest,anim={},unutra,brisi},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
SymDIRECT::maxiter="Dosegnut maksimalan broj iteracija `1`.";
n=Length[dom];
(* skaliranje funkcije na hiperkocku duljine stranice 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* početni centar i promjeri *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
unutra={test[c[[1]],d[[1]]]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};Nf=Nf+1;
TMin={c[[1]],fc[[1]]};
(* GLAVNA PETLJA *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* KORAK 1: traženje potencijalno optimalnih hiperkvadara *)
dmaxrazni=Union[dmax];
(* svi podaci grupirani po d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* podaci s najmanjim f grupirani po d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* provjera ostalih uvjeta *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* KORAK 2: dijeljenje potencijalno optimalnih *)
brisi={};
Do[
(* smjerovi po kojima se dijeli *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];Nf=Nf+Length[novic];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
(*uklanjanje onih koji nisu unutra*)
If[unutra[[pot]]==0,
(*centralni pravokutnik može izaći van ili postati sasvim unutra*)
If[test[c[[pot]],d[[pot]]]==-1,AppendTo[brisi,{pot}],unutra[[pot]]=test[c[[pot]],d[[pot]]]];
temp={};noviunutra={};
Do[
temp2=test[novic[[i]],novid[[i]]];
If[temp2==-1,AppendTo[temp,{i}],AppendTo[noviunutra,temp2]];
,{i,Length[novic]}];
novic=Delete[novic,temp];
novid=Delete[novid,temp];
novifc=Delete[novifc,temp];
,
(*ako je pravokutnik koji se dijeli sav unutra, kad se podjeli, svi manji pravokutnici će također biti unutra*)
noviunutra=Table[1,{Length[novic]}];
];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
unutra=Join[unutra,noviunutra];
(*update trenutnog minimuma*)
temp=Position[novifc,Min[novifc]][[1,1]];
If[novifc[[temp]]<TMin[[2]],TMin=N[{novic[[temp]],novifc[[temp]]},10]];
,{pot,pot2}];
(*brisanje onih centralnih koji su izašli van*)
c=Delete[c,brisi];
d=Delete[d,brisi];
dmax=Delete[dmax,brisi];
fc=Delete[fc,brisi];
unutra=Delete[unutra,brisi];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[SymDIRECT::maxiter,OptionValue[MaxIter]];];
(* slika za 2D slučaj *)
(*If[n==2 && iter<100,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Tooltip[Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10],unutra[[i]]]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}];
(* vraćanje u početnu domenu *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* ispis slika za 2D slučaj *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Data set:  5 Clusters

(* Data construction *)
SeedRandom[13]
n=2;q=2;k=5; nr=200;
boja={Blue,Green,Orange,Brown,Purple,Magenta}
prag=0.5;
p=.7;  kappa=Quantile[ChiSquareDistribution[n],p];
Q[th_]:={{Cos[th],-Sin[th]},{Sin[th],Cos[th]}}
centers={{2.5,2},{4.,5.},{4,7.5},{6.5,5.5},{8,2}};
cov={{{0.7,.1},{.1,0.2}},{{2.0,.01},{.01,.3}},{{1.5,-.2},{-.2,0.1}},{{0.4,.7},{.7,2}},{{0.3,-0.2},{-0.2,0.3}}};
Mk={150,250,140,180,130};Print["Total = ",Total[Mk]];
Do[
Print["C=",centers[[j]],"; cov=",cov[[j]],"; mj=",Mk[[j]]]
,{j,k}]
data=Table[RandomVariate[MultinormalDistribution[centers[[j]],cov[[j]]],Mk[[j]]],{j,k}];
P0=Table[data[[j]],{j,k}];
slcenters=ListPlot[centers,PlotStyle->{Brown,PointSize[.04],Opacity[.8]}];
A=Flatten[data,1];
AG=Table[{Floor[Min[A[[All,s]]]],Ceiling[Max[A[[All,s]]]]},{s,n}]
AG={{0,10},{0,10}};
slA=sl0=ListPlot[A,PlotStyle->{PointSize[.02]},AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
m=Length[A]; 
		(* Membership Matrix *)
mk=Table[Sum[Mk[[s]],{s,j}],{j,k}];mk=Prepend[mk,0];
U=Table[0,{i,m},{j,k}];
Do[
Do[
If[mk[[j]]+1≤ i≤  mk[[j+1]],U[[i,j]] =1]
,{j,k}]
,{i,m}]
Print["F = ",MFz[A,centers,U,cov]]
sldata=ListPlot[data,PlotStyle->{PointSize[.02]}];
If[n==2,slA1=Show[sldata,slcenters,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,Frame->True,ImageSize->Small]]
z0={{1,1},{4,4},{5,7},{7,7},{9,1}};
slz0=ListPlot[z0,PlotStyle->{Gray,Opacity[.7],PointSize[.04]}];
sl2=Show[slA,AspectRatio->Automatic,Frame->True,GridLines->Automatic,PlotRange->AG,ImageSize->Small]

MFz[A,centers,U_,cov]
MF[A,centers,U_,cov]

Data set:  3 Clusters

SeedRandom[13]
n=2;q=2.75;k=3; np=20;nr=2;sigma=0.1; tab=Table[i,{i,10}];
centers={{3.5,5.},{4,7.5},{7.5,5.5}};
cov={{{1.25,.05},{.05,.3}},{{1.3,-.3},{-.3,0.2}},{{0.4,.7},{.7,2}}};
AG={{0,10},{0,10}};pod=pi=Table[{},{j,k}];  Nr=Floor[np Area[Rectangle[{AG[[1,1]],AG[[2,1]]},{AG[[1,2]],AG[[2,2]]}]]];
data=RandomVariate[UniformDistribution[AG],Nr];
slcenters=ListPlot[centers,PlotStyle->{Brown,PointSize[.04],Opacity[.8]}];
Do[
Do[
If[dM[data[[i]],centers[[j]],5cov[[j]]]≤ 1,pod[[j]]=AppendTo[pod[[j]],data[[i]]]
]
,{i,Nr}]
,{j,k}]
Fig=ListPlot[pod,PlotRange->AG,AspectRatio->Automatic,Frame->True,GridLines->{tab,tab},ImageSize->Small]
Do[
Do[
pi[[j]]=AppendTo[pi[[j]],RandomVariate[MultinormalDistribution[pod[[j,s]],sigma IdentityMatrix[2]],nr]]
,{s,Length[pod[[j]]]}]
,{j,k}]
Pod=Table[Union[pod[[j]],Flatten[pi[[j]],1]],{j,k}];
sl0=ListPlot[Pod,PlotRange->AG,AspectRatio->Automatic,Frame->True,GridLines->Automatic,ImageSize->Small];
Show[sl0,slcenters]
A=Flatten[Pod,1]; m=Length[A];
slA=ListPlot[A,PlotRange->AG,AspectRatio->Automatic,Frame->True,GridLines->Automatic,ImageSize->Small]
Mk=Table[Length[Pod[[j]]],{j,k}];
		(* Membership Matrix *)
mk=Table[Sum[Mk[[s]],{s,j}],{j,k}];mk=Prepend[mk,0];U=Table[0,{i,m},{j,k}];
Do[
Do[
If[mk[[j]]+1≤ i≤  mk[[j+1]],U[[i,j]] =1]
,{j,k}]
,{i,m}]
Print["F = ",MFz[A,centers,U,cov]]
Print["Veličine klastera: ",Mk,";  |A| = ",m]

Call GK c - mens Data set:  5 Clusters

time=Timing[minGK=GKcmeans[A,z0,1]];
			(* Results and Figures *)
Print["Time=",time[[1]],"; It-",minGK[[4]],": centers=",minGK[[1]], "; FN = ",MFz[A,minGK[[1]],minGK[[2]],minGK[[3]]],"; Hdist=",Hdist[centers,minGK[[1]]]];
(* Rand Indexes *)
Print["Time - Rand and Jaccard (Frigue2007): ",RandFrigue[U,minGK[[2]]]//Timing]
Print["Time - Hullermeier: ",Hullermeier[U,minGK[[2]]]//Timing]
Print["Contingence matrix (dim): ",dim]
(* Attention: generally the matrix CC is of the tipe (k,l) *)
CC=Transpose[U].minGK[[2]]; dim=Dimensions[CC];
CC2=Table[Append[CC[[i]],Total[CC[[i]]]],{i,dim[[1]]}]; (* added a column *)
CC2=Append[CC2,Table[Sum[CC2[[i,j]],{i,dim[[1]]}],{j,dim[[2]]+1}]]; (* added a row *)
CCN=CC2;
Do[
Do[
CCN[[i,j]]=100CC2[[i,j]]/CC2[[i,dim[[2]]+1]]
,{j,dim[[2]]+1}]
,{i,dim[[1]]}]
Do[CCN[[dim[[1]]+1,j]]=Sum[CCN[[i,j]],{i,dim[[1]]}],{j,dim[[2]]}]; 
Print["C=",MatrixForm[Round[CC2,-.1]],"; C(%)=",MatrixForm[Round[CCN,-.1]]]
(* Figures *)
If[n==2,
slcenN=ListPlot[minGK[[1]],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.04],Opacity[.6]}];
cl=Table[minGK[[2,All,j]],{j,k}];
boje={1,3,6,8,1.5}/10; 
(* Nijanse *)
 nijansa={3,3,3,3,1};
slrez=Show[Table[Table[ListPlot[{A[[i]]},PlotStyle->{Hue[boje[[j]]],Opacity[cl[[j,i]]/nijansa[[j]]]}],{i,m}],{j,k}],slcenters,slcenN,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[slrez];
GraphicsGrid[{{slA1,slrez}},ImageSize->Medium]
]

Samo završna slika - dugo traje

slcenN=ListPlot[minGK[[1]],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.04],Opacity[.6]}];
cl=Table[minGK[[2,All,j]],{j,k}];
boje={1,3,6,8,1.5}/10; 
(* Nijanse *)
 nijansa={3,3,3,3,1};
slrez=Show[Table[Table[ListPlot[{A[[i]]},PlotStyle->{Hue[boje[[j]]],Opacity[cl[[j,i]]/nijansa[[j]]]}],{i,m}],{j,k}],slcenN,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[slrez];
GraphicsGrid[{{slA1,slrez}},ImageSize->Medium]

Call GK c - mens Data set:  3 Clusters

z0={{3,5},{4,7},{7,6}};
tab=Table[i,{i,10}];
time=Timing[minGK=GKcmeans[A,z0,1]];
			(* Results and Figures *)
Print["Time=",time[[1]],"; It-",minGK[[4]],": centers=",minGK[[1]], "; FN = ",MFz[A,minGK[[1]],minGK[[2]],minGK[[3]]],"; Hdist=",Hdist[centers,minGK[[1]]]];
(* Rand Indexes *)
Print["Time - Rand and Jaccard (Frigue2007): ",RandFrigue[U,minGK[[2]]]//Timing]
Print["Time - Hullermeier: ",Hullermeier[U,minGK[[2]]]//Timing]
(* Attention: generally the matrix CC is of the tipe (k,l) *)
CC=Transpose[U].minGK[[2]]; dim=Dimensions[CC];
Print["Contingence matrix (dim): ",dim]
CC2=Table[Append[CC[[i]],Total[CC[[i]]]],{i,dim[[1]]}]; (* added a column *)
CC2=Append[CC2,Table[Sum[CC2[[i,j]],{i,dim[[1]]}],{j,dim[[2]]+1}]]; (* added a row *)
CCN=CC2;
Do[
Do[
CCN[[i,j]]=100CC2[[i,j]]/CC2[[i,dim[[2]]+1]]
,{j,dim[[2]]+1}]
,{i,dim[[1]]}]
Do[CCN[[dim[[1]]+1,j]]=Sum[CCN[[i,j]],{i,dim[[1]]}],{j,dim[[2]]}]; 
Print["C=",MatrixForm[Round[CC2,-.1]],"; C(%)=",MatrixForm[Round[CCN,-.1]]]
(* Figures *)
If[n==2,
slcenN=ListPlot[minGK[[1]],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.04],Opacity[.6]}];
slrez=Show[slA,slcenters,slcenN,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[slrez];
]

Razdvajanje elemenata klastera

tab=Table[i,{i,10}];
Do[
yy=minGK[[2,All,j]];
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real}},Evaluate[Sum[Min[{(c1-yy[[i]])^2,(c2-yy[[i]])^2,(c3-yy[[i]])^2}],{i,Length[yy]}]]];
Nf=0;
dom=Table[{0,1},{3}];
time=Timing[
min=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->80,Eps->10^(-4)]
];
cen=Sort[min[[1]]];
a1=Mean[{cen[[1]],cen[[2]]}];
a2=Mean[{cen[[2]],cen[[3]]}];
Print["Centri: ",cen,"; Nf2=",Nf];
Print["ostatak: ",{0,a1},"; border: ",{a1,a2},"; core: ",{a2,1}];
Data=Table[{},{i,k}];
Do[
If[yy[[i]]<a1,Data[[1]]=Append[Data[[1]],A[[i]]],
If[a1≤ yy[[i]]<a2,Data[[2]]=Append[Data[[2]],A[[i]]],
Data[[3]]=Append[Data[[3]],A[[i]]]]]
,{i,Length[yy]}];
Print["|Data|:",DD=Table[Length[Data[[i]]],{i,k}],"; Total = ",Total[DD]];
cluster=Data[[3]];slcl=ListPlot[cluster,PlotStyle->{Brown,Opacity[1]}];
border=Data[[2]];slbor=ListPlot[border,PlotStyle->{Orange,Opacity[.5]}];
sl1=Show[slcl,slbor,slcenN,AspectRatio->Automatic,GridLines->{tab,tab},PlotRange->AG,AxesOrigin->{0,0},ImageSize->Small];
(*  *)
Dist[T_,Circ_]:=Module[{},Indikator=1;Return[Abs[Sqrt[dM[T,Circ[[1]],Circ[[3]]]]-Circ[[2]]]^2]
                                             ];
(* LAD pristup *)
C0=Median[border];  mb=Length[border];
cov=Sum[KroneckerProduct[(border[[s]]-C0),(border[[s]]-C0)],{s,mb}]/mb;
r0=Median[Table[Sqrt[dM[C0,border[[s]],cov]],{s,mb}]];   (* LADd *)
(* LAD *)
Print["C0=",C0];
F[pp_,qq_,rr_,al1_,al2_,al3_]:=Sum[Dist[border[[s]],{{pp,qq},rr,{{al1,al2},{al2,al3}} }],{s,mb}];
T=Timing[min=FindMinimum[{F[pp,qq,rr,al1,al2,al3],rr>0&&al1>0&&al3>0},{{pp,C0[[1]]},{qq,C0[[2]]},{rr,r0},{al1,cov[[1,1]]},{al2,cov[[1,2]]},{al3,cov[[2,2]]}}]];
Kr={{pp/.min[[2]],qq/.min[[2]]},rr/.min[[2]],{{al1/.min[[2]],al2/.min[[2]]},{al2/.min[[2]],al3/.min[[2]]}}};
slC1=ListPlot[{Kr[[1]]},PlotStyle->{Brown,PointSize[.04],Opacity[.4]}];
sle=ContourPlot[dM[Kr[[1]],{x,y},Kr[[3]]]==Kr[[2]]^2,
{x,0,10},{y,0,10},ContourStyle->{Red}];
sl2=Show[slbor,sle,slC1,slcenN,Axes->True,AspectRatio->Automatic,AxesOrigin->{0,0},PlotRange->{{0,10},{0,10}},GridLines->{tab,tab},ImageSize->Small];
Print[GraphicsGrid[{{Fig,sl1,sl2}},ImageSize->Large]]
,{j,k}]

Call Mahalanobis Fuzzy Incremental algorithm

mf=MFInc[A,{Mean[A]},5,0];
					(* Indexes and Graphs *)
Print["MFXB: ",Round[mf[[2,1]],-.005],"\nMFCH: ",Round[mf[[2,2]],-.005],"\nMFDB: ",Round[mf[[2,3]],-.005],"\nMFHV: ",Round[mf[[2,4]],-.005],"\nMFNew: ",Round[mf[[2,5]],-.000005]];Print[GraphicsGrid[{{ListPlot[Table[{i+1,mf[[2,1,i]]},{i,Length[mf[[2,1]]]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],
ListPlot[Table[{i+1,mf[[2,2,i]]},{i,Length[mf[[2,2]]]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],
ListPlot[Table[{i+1,mf[[2,3,i]]},{i,Length[mf[[2,3]]]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],
ListPlot[Table[{i+1,mf[[2,4,i]]},{i,Length[mf[[2,4]]]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic],
ListPlot[Table[{i+1,mf[[2,5,i]]},{i,Length[mf[[2,5]]]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic]}},ImageSize->Large]]
				(* Rand and Jaccard and Hausdorff *)
Print["Rand: ",Round[mf[[2,6]],-.005],"\nJaccard: ",Round[mf[[2,7]],-.005],"\nHausdorff: ",Round[mf[[2,8]],-.005]]
						(* Time *)
Print["Time(DIRECT) = ",mf[[3,1]],";  Time(GKc-means) = ",mf[[3,2]]]

Saving the data

SetDirectory[NotebookDirectory[]]; 
Export["PIP.txt",mf[[1]]]
Export["indexes.txt",mf[[2]]]

For the text

tab={2,3,4,5,6};
ListPlot[Table[{i+1,mf[[2,1,i]]},{i,Length[mf[[2,1]]]}],Joined->True,AxesOrigin->{1,0},Ticks->{Automatic,{.04,.08,.12,.16}},GridLines->{tab,{.04,.08,.12,0.16,.20}},ImageSize->Small]
ListPlot[Table[{i+1,mf[[2,2,i]]},{i,Length[mf[[2,2]]]}],Joined->True,AxesOrigin->{1,0},PlotRange->{0,2500},GridLines->{tab,{500,1000,1500,2000,2500}},ImageSize->Small]
ListPlot[Table[{i+1,mf[[2,3,i]]},{i,Length[mf[[2,3]]]}],Joined->True,AxesOrigin->{1,.1},Ticks->{Automatic,{.2,.3,.4}},GridLines->{tab,{.2,.3,.4}},PlotRange->{.1,.4},ImageSize->Small]
ListPlot[Table[{i+1,mf[[2,4,i]]},{i,Length[mf[[2,4]]]}],Joined->True,AxesOrigin->{1,1},GridLines->{Automatic,{1.5,2.5,3.5}},PlotRange->{1,4},ImageSize->Small]
ListPlot[Table[{i+1,mf[[2,5,i]]},{i,Length[mf[[2,5]]]}],Joined->True,AxesOrigin->{0,0},PlotRange->{0,0.01},GridLines->{Automatic,{0,0.002,0.004,0.006,0.008,0.01}},ImageSize->Small]

Optimal partition

z0={{2.5612312296761885`,2.000265017358591`},{7.87789950000471`,2.123341045967437`},{6.400295638968773`,5.628445590522856`},{3.839609529981384`,7.513944963062942`},{3.3594818219995957`,4.921384360910297`}};
time=Timing[minGK=GKcmeans[A,z0,0]];
Print["cen: ",minGK[[1]],"\ncov: ",minGK[[3]],"\nMFz = ",MFz[A,minGK[[1]],minGK[[2]],minGK[[3]]],": MF = ",MF[A,minGK[[1]],minGK[[2]],minGK[[3]]]];slcenN=ListPlot[minGK[[1]],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.04],Opacity[.6]}];
cl=Table[minGK[[2,All,j]],{j,k}];
boje={1,3,6,8,1.5}/10; 
(* Nijanse *)
 nijansa={3,3,3,3,1};
slrez=Show[Table[Table[ListPlot[{A[[i]]},PlotStyle->{Hue[boje[[j]]],Opacity[cl[[j,i]]/nijansa[[j]]]}],{i,m}],{j,k}],slcenters,slcenN,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[slrez];
GraphicsGrid[{{slA1,slrez}},ImageSize->Medium]
Print["Contingence matrix (dim): ",dim]
(* Attention: generally the matrix CC is of the tipe (k,l) *)
CC=Transpose[U].minGK[[2]]; dim=Dimensions[CC];
CC2=Table[Append[CC[[i]],Total[CC[[i]]]],{i,dim[[1]]}]; (* added a column *)
CC2=Append[CC2,Table[Sum[CC2[[i,j]],{i,dim[[1]]}],{j,dim[[2]]+1}]]; (* added a row *)
CCN=CC2;
Do[
Do[
CCN[[i,j]]=100CC2[[i,j]]/CC2[[i,dim[[2]]+1]]
,{j,dim[[2]]+1}]
,{i,dim[[1]]}]
Do[CCN[[dim[[1]]+1,j]]=Sum[CCN[[i,j]],{i,dim[[1]]}],{j,dim[[2]]}]; 
Print["C=",MatrixForm[Round[CC2,-.1]],"; C(%)=",MatrixForm[Round[CCN,-.1]]]