Mahalanobis weighted Incremental algorithm

Once the set of optimal partitions is obtained, one should activate the M-Indexes[] module to get the appropriate Mahalanobis indexes
    
    Input: data set 𝒜 ⊂ (ℝ^n)data weights w initial center(s) cen. number of iterations K. Each additional center is obtained solving a GOP using Compile.In case when n>2, one has to adapt the minimizing function.
Output: set of partitions each partition consists of clusters with corresponding weights associated to their elements).Graphical depictions of each iteration are presented only for n=2;;
          Submodules: dM[], MWF[], MFmin[], Hdist[], DIRECT[], MWKMeans[]

(* Ako je poznata originalna particija PI0, možemo je usporediti s particijom dobivenom u brk-koraku *)
MInc[w_,A_,cen_,KK_,Ind_:0]:=Module[{centri,centriN,m=Length[A],n=Length[A[[1]]],k,dom,Scov,min,PI1},
Print["n=",n];
(* Initialization *)
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
Scov=Table[IdentityMatrix[n],{i,Length[cen]}]; 
sol=MWKMeans[w,A,cen]; 
(*cen=sol[[2]];*)
centri=sol[[2]];
brk0=Length[centri];
PART={};
					(* Iterativni postupak *)
Do[
If[Ind≠0,Print["\nNa postojecih ",Length[centri]," centara dodaje se 1 novi"]];
k=Length[centri];
						(* Compile 
FCompile=Compile[{{u,_Real},{v,_Real}},
Evaluate[Sum[w[[i]]Min[Append[Table[Norm[A[[i]]-centri[[s]]]^2,{s,Length[centri]}],Norm[A[[i]]-{u,v}]^2]],{i,m}] ]];	
F2[u_,v_]:=Apply[FCompile,{u,v}];
*)
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[A[[i]]-centri[[s]]]^2,{s,Length[centri]}],Norm[A[[i]]-{u,v}]^2]],{i,m}];
dom={xG,yG};
	(* DIRECT *)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];If[Ind≠0,Print["DIRECT: ",min]];
centriN=Append[centri,min[[1]]];
If[Ind≠0,Print["Novi optimirani centar: ",centriN]];
(* Ovdje se vrijednost funkcije cilja ne moze racunati kao MWFmin jer jos nemam Scov[[k+1]] *)
k=Length[centriN];
					(* Poziv Mahalanobis k-means *)
t1=Timing[sol=MWKMeans[w,A,centriN]];
If[Ind≠0,Print["Novi MWKMeans profinjeni centar: ",sol[[2]]]];
Scov=sol[[3]];
If[Ind≠0,Print["MWKMeans iteracija = ",sol[[4]],"; MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]]];
(* Ako je poznata originalna particija, ovdje ju možemo usporediti s particijom dobivenom u brk-koraku:
PI1=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
If[brk≤Length[PI0],rand=Rand[PI1,PI0],rand=Rand[PI0,PI1]];
Print["Confusion Matrix = ",rand[[1]]//MatrixForm, "; Rand = ",rand[[2,1]],"; Jaccard = ",rand[[2,2]]];
*)
If[Ind≠0,Print["Time (DIRECT+WKM) = ",minp[[1]],"+",t1[[1]]," = ",minp[[1]]+t1[[1]]]];
Print[brk,"-particija"];
     (* FIGUREE *)
k=Length[sol[[1]]];
If[n==2,
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.05],Opacity[.7]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{Table[PointSize[.015],{j,k}]},PlotStyle->{Table[Hue[j/k],{j,k}]}]
];
							(* Elipse *)
If[n==2,
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}]
,{j,k}];
sl=Show[Table[slcon[j],{j,k}],ImageSize->Small];
sl2=Print[Show[slPI,slc,sl,GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]];	
];		 
centri=sol[[2]];     
PART=Append[PART,sol[[1]]]; 
,{brk,KK}];
PART
]