MIncAdd - clustering

Input: data set B⊂ (ℝ^3)x (Z^2)centers (Subscript[cen, 1]),(Subscript[cen, 2]),(Subscript[cen, 3]) ∈ (([0,1])^3), weights (Subscript[w, 1]),⃛,(Subscript[w, m]) >0,number of iterations KK.Each additional center is obtained solving a GOP
using Compile;
Output: Optimal partition,  centers, covariance matrices and eigenvalues and eigenvectors of p-clusters.
Submodules: dM[],MWF[], MWKMeans2[], MWKMeansAdd[], IASymDIR[], ASymDIR[], DIRECT[], Hdist[]

MIncAdd[w_,B_,cen_,KK_,Ind_:0]:=Module[{centri,dom,Scov,min},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,Length[B]}]; 
AG=Table[{0,1},{i,3}];
				(* Initialization *)
centri=cen;   
Do[
If[brk>1,
If[Ind≠0,Print[Style["Iteracija ",Red,Bold,12],Style[brk,Red,Bold,14]," Na postojećih ",Length[centri]," centara dodaje se 1 novi"]];
k=Length[centri];  
					(* DIRECT Compile *)
time1=Timing[
Clear[u1,u2,u3,F2];
FCompile=Compile[{{u1,_Real},{u2,_Real},{u3,_Real}},
Evaluate[Sum[w[[i]]Min[Append[Table[EuclideanDistance[Pod[[i,2]],centri[[s]]],{s,Length[centri]}],EuclideanDistance[Pod[[i,2]],{u1,u2,u3}]]],{i,Length[B]}] ]];	
F2[u1_,u2_,u3_]:=Apply[FCompile,{u1,u2,u3}];	
dom=AG;					
min=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)];
];
centriN=Append[centri,min[[1]]],
centriN=cen];
							(*  *)
If[Ind≠0,Print["Novi optimirani centri: ",centriN]];
k=Length[centriN];
			(* Poziv Mahalanobis k-means *)
time2=Timing[ 
If[k==2,sol=MWKMeans2[w,B,centriN,Ind],
 sol=MWKMeansAdd[w,B,centriN,Ind]]
(* sol=MWKMeansAdd[w,B,centriN,Ind]  *)
                    ]; 
tab=Table[{sol[[1,j]],sol[[2,j]],sol[[3,j]],Norm[sol[[2,j]]],j},{j,Length[sol[[2]]]}];
sort=Sort[tab,#1[[4]]<#2[[4]]&];
sol[[1]]=sort[[All,1]];sol[[2]]=sort[[All,2]];sol[[3]]=sort[[All,3]];
If[Ind≠0,Print["Novi MWKMeans profinjeni centri: ",sol[[2]]]];
Scov=sol[[3]];
If[Ind≠0,Print["MWKMeans iteracija: iter=",sol[[4]],"; MWF = ",FF=MWF[sol[[1]],sol[[2]],sol[[3]]],"; |pi|: ",mjx=Table[Length[sol[[1,j]]],{j,Length[sol[[1]]]}]]];
                                                             (* Construction of MWKM image *)
If[Ind2≠0,
dol=sol;
Do[dol[[1,2,s,2]]={1,1,1},{s,Length[dol[[1,2]]]}];PIP=Flatten[Table[Table[{dol[[1,j,s,2]],dol[[1,j,s,3]]},{s,Length[dol[[1,j]]]}],{j,Length[dol[[2]]]}],1];
sl=SortBy[PIP,Last];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print["MWKM",Length[sol[[2]]],"-image.",  " Optimalni centri: ",sol[[2]],"; F(MWF) = ",MWF[sol[[1]],sol[[2]],sol[[3]]],"\np-clusters size: ",{Length[sol[[1,1]]],Length[sol[[1,2]]]},"; iter=",sol[[4]]];
Print[slN=Image[slikaN,ImageSize->Small]];
        ];
                                                           (* pPI0, cov, Eigenvalues *)
If[Ind1≠0,
time3=Timing[
pPI0=Table[{},{j,Length[sol[[1]]]}];
PI0=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];  (* Clusters bez koordinata pixela *)
cc=Table[Mean[PI0[[j]]],{j,Length[sol[[1]]]}];
cov=Table[Covariance[PI0[[j]]],{j,Length[sol[[1]]]}];
(*ell=Table[Ellipsoid[cc[[j]], 11.3449  cov[[j]]],{j,Length[PI0]}];*)(* p=.99 *)
pell=Table[Ellipsoid[cc[[j]], kappa  cov[[j]]],{j,Length[sol[[1]]]}];
Do[
Do[
If[PI0[[j,s]]∈pell[[j]] ,pPI0[[j]]=Append[pPI0[[j]],PI0[[j,s]]]  ];
,{s,Length[PI0[[j]]]}]
,{j,Length[sol[[1]]]}];
cov=Table[Covariance[pPI0[[j]]],{j,Length[sol[[1]]]}];
ei=Table[Eigenvalues[cov[[j]]],{j,Length[sol[[1]]]}];
                                                    (* Ellipsoids *)
If[Ind3!=0,
slell=Show[Table[Graphics3D[{RGBColor[cc[[j]]],Opacity[0.4],pell[[j]]}],{j,Length[sol[[2]]]}],Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"}];
sltoc=Table[ListPointPlot3D[pPI0[[j]],PlotStyle->{RGBColor[cc[[j]]],PointSize[.008],Opacity[.8]}],{j,Length[cc]}];
Print[Show[slell,sltoc, ImageSize->150]];
  ];
];
];
					(* Ellipsoids + Eigenvalues *)
If[Ind1==0,
time3=Timing[
PI0=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];                  (* Clusters bez koordinata pixela *)                                    
cc=Table[Mean[PI0[[j]]],{j,Length[sol[[1]]]}];
cov=Table[Covariance[PI0[[j]]],{j,Length[sol[[1]]]}];
ei=Table[Eigenvalues[cov[[j]]],{j,Length[sol[[1]]]}];
ell=Table[Ellipsoid[cc[[j]], 11.3449  cov[[j]]],{j,Length[sol[[1]]]}];
pell=Table[Ellipsoid[cc[[j]], kappa  cov[[j]]],{j,Length[sol[[1]]]}];
                                              (* Ellipsoids *)
If[Ind3!= 0,slell=Show[Table[Graphics3D[{RGBColor[cc[[j]]],Opacity[0.4],pell[[j]]}],{j,Length[sol[[1]]]}],Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"}];
sltoc=Table[ListPointPlot3D[PI0[[j]],PlotStyle->{RGBColor[cc[[j]]],PointSize[.008],Opacity[.8]}],{j,Length[cc]}];
Print[Show[slell,sltoc, ImageSize->150]];
   ];
];
];	
centri=sol[[2]];
If[Ind5≠0,
Print[brk+1,Style["-partition of the set data2:",13,Blue]];
Print["Centroids: ",centri, "; MWF= ",sol[[5]]];
Print[Style["Eigenvalues: ",13,Blue],ei];
EIG[[brk]]=Append[EIG[[brk]],ei];
];
,{brk,KK}];
{sol[[1]],sol[[2]],sol[[3]],ei}]