Recognizing of a linearly non - separable datasets
1-3-2023

Modules

Data construction

(* Data coming from a Circle *)
DataCircle[cen_,r_,m_,sigma_,al_,be_]:=Module[{A={},t,T},
Do[t=RandomReal[{al,be}];
T=cen+r {Cos[t],Sin[t]};
A=Append[A,RandomReal[MultinormalDistribution[T,sigma IdentityMatrix[2]],{1}][[1]]];,{ii,m}];
A
]

(* Data coming from a Line *)
DataLine[cA_,cB_,broj_,sigma_]:=Module[{Mk=broj Norm[cA-cB],pod={},lam,T},
Do[
lam=s/Mk; 
T=lam cA+(1-lam)cB;
pod=Append[pod,RandomReal[MultinormalDistribution[T,sigma IdentityMatrix[2]],{1}][[1]]];
,{s,Mk}];
pod
]

(* Data coming from a Ellipse with center OO={0,0} *)
DataEll[a_,b_,m_,sigma_]:=Module[{tt,CC,distr,tocke1,},
CC=1/(a*Integrate[Sqrt[1-(1-b^2/a^2)*Cos[x]^2],{x,0,2Pi}])//N;
distr=ProbabilityDistribution[CC*a*Sqrt[1-(1-b^2/a^2)*Cos[x]^2],{x,0,2Pi}];
tt=RandomVariate[distr,m];
tocke1=Transpose[{a*Cos[tt],b*Sin[tt]}];
tockeoblak=tocke1+RandomVariate[NormalDistribution[0,sigma],{m,2}]
]

Auxiliary programs

Weighted k - means   for the data with additional information with Break

WKMeansDop[w_,B_,z_,Ind_:0]:=Module[{m=Length[B],k=Length[z],n=Length[B[[1]]],c,c0,Pod,imin,PI,tab,It=0},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,m}]; 
c0=z; 
PI=Table[{},{i,k}]; 
			(* Minimal distance principle *)
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
(*Print["Initial centers: ",z,"; F = ",WF[PI]];*)
						(* New centers *)
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition
				in the case n=2 - can be extended to an arbitrary n *)
If[Length[c]<Length[c0],
Print["There was a case of reducing the number of clusters in the partition"];Break[] ];
Chop[Norm[c-c0]]≠0&&It<25,c0=c;It=It+1;
				(* New partition *)
PI=Table[{},{i,k}]; 
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
               ];
{PI,c0,WF[PI],It}]

Weighted objective functions

Input: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} of the set 𝒜 ⊂ (ℝ^n), (in case of WFmin[] only the set 𝒜), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n);
Output: weighted objective function values: ℱ(Π), \[FormalCapitalF](c), 𝒢(Π);
Submodules: none needed

WF[PI_]:=Module[{m=Length[Flatten[PI,1]],c},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Return[Sum[PI[[j,s,1]]*Norm[PI[[j,s,2]]-c[[j]]]^2,{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

WFmin[PI_]:=Module[{A=Flatten[PI,1][[All,2]],k=Length[PI]},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Sum[Flatten[PI,1][[All,1]][[i]]Min[Table[Norm[c[[j]]-A[[i]]]^2,{j,k}]],{i,Length[A]}]
]

WG[PI_]:=Module[{m=Length[Flatten[PI,1]],c,set,cc},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
set=Flatten[PI,1];
cc=Mean[WeightedData[set[[All,2]],set[[All,1]]]];
Return[Sum[Length[PI[[j]]]*Norm[cc-c[[j]]]^2,{j,Length[PI]} ]]
];

LinConstrDIRECTe

Modul for solving the global optimization problem for a Lipschitz continuous function with linear constraints;
Input: ordinary or compiled function, Domain , {Eps->10^(-4),MaxIter->20,DMin->10.^(-2)};
Region ℛ1 given by  A1x≤ b1, Initial points pts;
Output: optimal variables values, optimal objective function value , CPU-time;
Submodules: Test3[], Test3plus[]

(* Vrijednost aktualnog minimuma određuje se pomoću vrijednosti funkcije u točki koja se nalazi na presjeku pravokutnika i područja R1, a koja se dobiva prema Strategy~1 ili Strategy~2;
Najprije se početna vrijednost TMin definira kao TMin={pts,fc[[1]]}; gdje je pts izabrana točka u R1; 
Prametar iter je izvađen iz popisa internih varijabli;
 *)
Options[LinConstrDIRECTe]={Eps->0,MaxIter->50,DMin->10.^(-3)};
LinConstrDIRECTe[f_,dom_,OptionsPattern[]]:=
Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,temp2,ind,novic,novid,novifc,noviunutra,red,TMin,dtest,anim={},unutra,brisi,sign,rec,pom},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
LinConstrDIRECTe::maxiter="Dosegnut maksimalan broj iteracija `1`.";
n=Length[dom];
iter=0;
tn={};
sign=Tuples[{-1,1},n];
(* scaling the function on the hypercube with length of the side 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* initial center and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
unutra={Test3[c[[1]],d[[1]]]};  
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@pts};Nf=Nf+1;
TMin={pts,fc[[1]]};   (* Print[iter,". ",TMin]; *)
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: searching for potentially optimal hyperrectangles *)
dmaxrazni=Union[dmax];
(* All data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* Data with smallest f grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* Checking other conditions *)
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];
(* STEP 2: dividing potentially optimal hyperrectangles *)
brisi={};
Do[
(* Dividing directions *)
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]}];
(* Eliminating those that are not inside *)
If[unutra[[pot]]==0,
(* Central hyperrectangle can go out or stay completely inside *)
pom=Test3plus[c[[pot]],d[[pot]],ff,TMin];
If[pom==-1,AppendTo[brisi,{pot}],unutra[[pot]]=pom;TMin=tmin];
temp={};noviunutra={};
Do[
temp2=Test3plus[novic[[i]],novid[[i]],ff,TMin];TMin=tmin;
If[temp2==-1,AppendTo[temp,{i}],AppendTo[noviunutra,temp2]];
,{i,Length[novic]}];
novic=Delete[novic,temp];
novid=Delete[novid,temp];
novifc=Delete[novifc,temp];,
(* If the divided hyperrectangle lies completely inside, all subhyperrectangles remain inside too *)
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 of actual minimum *)
If[novifc≠ {},temp=Position[novifc,Min[novifc]][[1,1]]]; 
(*Print["novifc: ",novifc,"\nMin: ",Min[novifc]];*)
(*
rec=Table[Table[novic[[temp,j]]+sign[[i,j]]*novid[[temp,j]]/2,{j,n}],{i,Length[sign]}]; 
Do[
If[RegionMember[ℛ1,rec[[j]]]==True&&ff@@rec[[j]]<TMin[[2]],TMin=N[{rec[[j]],ff@@rec[[j]]},10]],{j,Length[rec]}];
Mat=({{Join[A1,IdentityMatrix[n1],-IdentityMatrix[n1],IdentityMatrix[n1]], IdentityMatrix[m1+3n1]}})//ArrayFlatten;
vec=Join[b1,Table[1,{i,n1}],-novic[[temp]]+novid[[temp]]/2,novic[[temp]]+novid[[temp]]/2];
Table[If[vec[[i]]<0,vec[[i]]=-vec[[i]];Mat[[i]]=-Mat[[i]]],{i,Length[vec]}];
Mat=({{Mat, IdentityMatrix[m1+3n1]}})//ArrayFlatten;
lp=LinearProgramming[Join[Table[0,{i,4n1+m1}],Table[1,{i,m1+3n1}]],Mat,Table[{vec[[i]],0},{i,Length[vec]}]];
If[Norm[Take[lp,-(m1+3n1)]]==0,pt=Take[lp,n1];ft=Apply[ff,pt];Nf=Nf+1;,pt={};ft=+Infinity];
If[ft<TMin[[2]],TMin=N[{pt,ft}]] ; 
*)
(*tn=Append[tn,N[{novic[[temp]],novid[[temp]],novifc[[temp]]},5]];*)
,{pot,pot2}];  (* Print[iter,". ",TMin]; *)
(* Deletion of those hyperrectangles which went outside *)
c=Delete[c,brisi];
d=Delete[d,brisi];
dmax=Delete[dmax,brisi];
fc=Delete[fc,brisi];
unutra=Delete[unutra,brisi];
If[iter==OptionValue[MaxIter],Message[LinConstrDIRECTe::maxiter,OptionValue[MaxIter]];];
(* Figure in 2D case *)
(*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]}];
(* Back to the initial domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* Images printing in 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Test3: 
prvi modul Test3[] radi korištenjem vertices + LP, a drugi modul Test3[] radi korištenjem samo LP
Izmjenom redosljeda modula definira se aktualni modul

(* Test3 - LP *)
Test3[c_,dd_]:=Module[{n=Length[c],d=dd/2,tt=-1,sig,rec},
LP=LP+1;
strategy=1;
Mat=({{Join[A1,IdentityMatrix[n1],-IdentityMatrix[n1],IdentityMatrix[n1]], IdentityMatrix[m1+3n1]}})//ArrayFlatten;
(* Provjera je li presjek pravokutnika i R1 neprazan: koristi se LP *)
vec=Join[b1,Table[1,{i,n1}],-c+d,c+d];
Table[If[vec[[i]]<0,vec[[i]]=-vec[[i]];Mat[[i]]=-Mat[[i]]],{i,Length[vec]}];
Mat=({{Mat, IdentityMatrix[m1+3n1]}})//ArrayFlatten;
If[Norm[Take[LinearProgramming[Join[Table[0,{i,4n1+m1}],Table[1,{i,m1+3n1}]],Mat,Table[{vec[[i]],0},{i,Length[vec]}]],-(m1+3n1)]]==0,tt=0 ; LP1=LP1+1 ;
    ];
(* Provjera je li cijeli pravokutnik u R1  u slučaju n≤3 *)
If[n≤3,
itt=0;
sig=Tuples[{-1,1},n];
rec=Table[Table[c[[j]]+sig[[k,j]]d[[j]],{j,n}],{k,Length[sig]}];
Do[
If[RegionMember[ℛ1,rec[[jj]]]==True,itt=itt+1]
,{jj,Length[rec]}];
If[itt==Length[rec],tt=1]; 
];
tt
]

Test3plus - LP

Uz funkciju Test3, dodatno se izračunava aktualnu vrijednost minimuma u točki pt koju daje LP;;
Input: center c and vector dd of half side-lengths of the  rectangle R[c,dd];
Output: -1 if R ⋂ℛ1 =Θ; 0 if R ⋂ℛ1 ≠ Θ; 1 if R ⊂ℛ1 
Submodules: none needed

(* 
U tu svrhu dodatno se poziva funkcija ff i aktualna vrijednost TMin;
Za dobivenu točku pt izračunava se ft=f@@pt; Ako je ft<TMin[[2]], update TMin;
Prije korištenja modula treba biti poznat region ℛ1, matrica A1 i vektor b1;
Na kraju ostaje provjera nalazi li se cijeli pravokutnik u ℛ1;
Ako se isključi provjera nalazi li se cijeli pravokutnik u ℛ1 CPU-time skrati se 10 puta (na G01: 13 varijabli); 
Ovo nema posljedice na traženje optimuma: to samo znači da će svi rectangles koji nastanu od ovakvog pravokutnika nadalje testirati. Ostane li ova promjena, ovi pravokutnici se ne bi trebali testirati.
*)
Test3plus[c_,dd_,ff_,TMin_]:=Module[{n=Length[c],d=dd/2,tt=-1,pt,ft,sig,rec},
LP=LP+1;
strategy=1;
Mat=({{IdentityMatrix[m1], ConstantArray[0,{m1,n1}], -A1, -IdentityMatrix[m1]}, {ConstantArray[0,{n1,m1}], IdentityMatrix[n1], IdentityMatrix[n1], ConstantArray[0,{n1,m1}]}});
Mat=ArrayFlatten[Mat];
vec=({{b1-A1.(c+d), 2 d}});
vec=Flatten[ArrayFlatten[vec][[1]]];
cc1=({{Table[0,{i,1,m1}], Table[0,{i,1,n1}], Table[0,{i,1,n1}] , Table[1,{i,1,m1}]}});
cc1=Flatten[ArrayFlatten[cc1][[1]]];
Table[If[vec[[i]]<0,vec[[i]]=-vec[[i]];Mat[[i]]=-Mat[[i]]],{i,Length[vec]}];
lp=LinearProgramming[cc1,N[Mat],N[Table[{vec[[i]],0},{i,Length[vec]}]]];
If[Norm[Take[lp,-m1]]==0,pt=c+d-Take[lp,{m1+n1+1,m1+2n1}];ft=ff@@pt;Nf=Nf+1;
If[ft<TMin[[2]],tmin=N[{pt,ft}],tmin=TMin];
LP1=LP1+1;
tt=0,  tmin=TMin
    ];
(* Provjera je li cijeli pravokutnik u R1 u slučaju n≤3 *)
If[n≤4,itt=0;
sig=Tuples[{-1,1},n];
rec=Table[Table[c[[j]]+sig[[k,j]]d[[j]],{j,n}],{k,Length[sig]}];
Do[
If[RegionMember[ℛ1,rec[[jj]]]==True,itt=itt+1]
,{jj,Length[rec]}];
If[itt==Length[rec],tt=1]; 
];
tt
]

Main program - Spectral Clustering Algorithms

S: a linearly non - separable dataset;
k: number of clusters;
Sigma: parameter [0.1, 5];
Alg: Spectar[1], SpectarLCDIR[2];
IT: number of iterations for SpectarLCDIR

SCAlg[S_,k_,Sigma_,Alg_,IT_]:=Module[{},
m=Length[S];
SDop=Table[{S[[i]],i},{i,m}]; 
A=Table[If[i==j,0.,Round[Exp[-Norm[SDop[[i,1]]-SDop[[j,1]]]^2/(2 Sigma^2)],10^(-6)]//N],{i,m},{j,m}];
DDD=Table[Total[A[[i]]],{i,m}];
 If[Product[DDD[[ii]],{ii,m}]≠0,  
L=Table[A[[i,j]]/Sqrt[DDD[[i]] DDD[[j]]]//N,{i,m},{j,m}];
eig=Eigensystem[L,k];
X=Transpose[eig[[2]]];
Y=Table[X[[i,j]]/Norm[X[[i]]],{i,m},{j,k}];
Y0=Table[{Y[[i]],i},{i,m}];
w=Table[1,{i,m}];
If[Alg==1,
				(* SPECTAR Algorithm *) 
(* z0 biram iz Y od međusobno što okomitijih redaka *)
time=Timing[z0=Table[0,{j,k}];
z0[[1]]=RandomChoice[Y];
Do[
Y1= Complement[Y,Table[z0[[j]],{j,i-1}]];
z0[[i]]=Y1[[1]]; sp0=Sum[Abs[z0[[i]].z0[[j]]],{j,i-1}];
Do[
x=Y1[[s]]; 
If[Sum[Abs[x.z0[[j]]],{j,i-1}]<sp0, z0[[i]]=x; sp0=Sum[Abs[z0[[i]].z0[[j]]],{j,i-1}]]
,{s,2,Length[Y1]}];
,{i,2,k}];
]; 
time1=Timing[z1=Orthogonalize[z0]];
				(* WKMeans *)
time2=Timing[sol=WKMeansDop[w,Y0,z1]];
sol0=sol[[1]];
z2=sol[[2]];
Print[Style["k = ",16,Blue],Style[Length[sol0],16,Blue],Style["; |S|=",16,Blue],Style[m,16,Blue],Style[";  F = ",16,Blue],Style[sol[[3]],16,Blue],Style["; Sigma = ",16,Blue],Style[Sigma,16,Blue],"\nThe highest eigenvalues: ",eig[[1]],"\nc = ",Round[sol[[2]],-.001], "; Iteracija: ",sol[[4]],"; CPU = ",time[[1]]];
Print["Ini. Appr.: ",Round[z0,-.0001]];
Print["Normalize: ",Round[z1,-.0001]];
Print["k-means: ",Round[sol[[2]],-.0001]];
Print["Time {Ini, Normalize, k-means}: ",{time[[1]],time1[[1]],time2[[1]]}," = ",time[[1]]+time1[[1]]+time2[[1]]];
				(* Konstrukcija particije rješenja *)
PI=Table[{},{j,k}];
Do[
ss=sol0[[j,All,3]];
PI[[j]]=Table[S[[s]],{s,ss}]
,{j,k}];
	];
If[Alg==2,
			(* z0 is determined with a few iteration of LCDIR[] *)
						(* Constraints *)
n=Length[Y[[1]]]; kn= k n;
A1=Table[0,{i, k-1},{j,kn}];
Do[
Do[A1[[i,j]]=-1; A1[[i,j+n]]=1;
,{j,n (i-1)+1,n i}]
,{i,k-1}];
b1=Table[0,{i,k-1}]; 
pts=Table[.25,{i,kn}];   (* Početna točka*)
Print["[A1;b1] = ",{MatrixForm[A1],MatrixForm[b1]}];
Print["{m1,n1} = ",{m1,n1}=Dimensions[A1],"; |Y| = ",Length[Y],"; pts = ",pts];
					(* Region with k*n variables *)
If[k==2,
ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4}];
FFCompile=Compile[{{x1,_Real},{x2,_Real},{x3,_Real},{x4,_Real}},Evaluate[Sum[Min[{Norm[{x1,x2}-Y[[i]]]^2,Norm[{x3,x4}-Y[[i]]]^2}],{i,m}]]]];
If[k==3,
ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4,x5,x6,x7,x8,x9}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4,x5,x6,x7,x8,x9}];
Print["R1: ",ℛ1];
FFCompile=Compile[{{x1,_Real},{x2,_Real},{x3,_Real},{x4,_Real},{x5,_Real},{x6,_Real},{x7,_Real},{x8,_Real},{x9,_Real}},Evaluate[Sum[Min[{Norm[{x1,x2,x3}-Y[[i]]]^2,Norm[{x4,x5,x6}-Y[[i]]]^2,Norm[{x7,x8,x9}-Y[[i]]]^2}],{i,m}]]]];
If[k==4,
ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16}];
FFCompile=Compile[{{x1,_Real},{x2,_Real},{x3,_Real},{x4,_Real},{x5,_Real},{x6,_Real},{x7,_Real},{x8,_Real},{x9,_Real},{x10,_Real},{x11,_Real},{x12,_Real},{x13,_Real},{x14,_Real},{x15,_Real},{x16,_Real}},Evaluate[Sum[Min[{Norm[{x1,x2,x3,x4}-Y[[i]]]^2,Norm[{x5,x6,x7,x8}-Y[[i]]]^2,Norm[{x9,x10,x11,x12}-Y[[i]]]^2,Norm[{x13,x14,x15,x16}-Y[[i]]]^2}],{i,m}]]]];
If[k==5,
ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25}];
FFCompile=Compile[{{x1,_Real},{x2,_Real},{x3,_Real},{x4,_Real},{x5,_Real},{x6,_Real},{x7,_Real},{x8,_Real},{x9,_Real},{x10,_Real},{x11,_Real},{x12,_Real},{x13,_Real},{x14,_Real},{x15,_Real},{x16,_Real},{x17,_Real},{x18,_Real},{x19,_Real},{x20,_Real},{x21,_Real},{x22,_Real},{x23,_Real},{x24,_Real},{x25,_Real}},Evaluate[Sum[Min[{Norm[{x1,x2,x3,x4,x5}-Y[[i]]]^2,Norm[{x6,x7,x8,x9,x10}-Y[[i]]]^2,Norm[{x11,x12,x13,x14,x15}-Y[[i]]]^2,Norm[{x16,x17,x18, x19,x20}-Y[[i]]]^2,Norm[{x21,x22,x23,x24,x25}-Y[[i]]]^2}],{i,m}]]]];
If[k==6,
ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36}];
FFCompile=Compile[{{x1,_Real},{x2,_Real},{x3,_Real},{x4,_Real},{x5,_Real},{x6,_Real},{x7,_Real},{x8,_Real},{x9,_Real},{x10,_Real},{x11,_Real},{x12,_Real},{x13,_Real},{x14,_Real},{x15,_Real},{x16,_Real},{x17,_Real},{x18,_Real},{x19,_Real},{x20,_Real},{x21,_Real},{x22,_Real},{x23,_Real},{x24,_Real},{x25,_Real},{x26,_Real},{x27,_Real},{x28,_Real},{x29,_Real},{x30,_Real},{x31,_Real},{x32,_Real},{x33,_Real},{x34,_Real},{x35,_Real},{x36,_Real}},Evaluate[Sum[Min[{Norm[{x1,x2,x3,x4,x5,x6}-Y[[i]]]^2,Norm[{x7,x8,x9,x10,x11,x12}-Y[[i]]]^2,Norm[{x13,x14,x15,x16,x17,x18}-Y[[i]]]^2,Norm[{x19,x20,x21,x22,x23,x24}-Y[[i]]]^2,Norm[{x25,x26,x27,x28,x29,x30}-Y[[i]]]^2,Norm[{x31,x32,x33,x34,x35,x36}-Y[[i]]]^2}],{i,m}]]]];
dom=Table[{-1,1},{kn}];
						(* Poziv LinConstrDIRECTe[] s IT iteracija *)
LP=LP1=Nf=0; 
time=Timing[min=LinConstrDIRECTe[FFCompile,dom,DMin->10.^(-3),MaxIter->IT, Eps->10^(-4)]];
If[RegionMember[ℛ1,min[[1]]]==True,Print["Dobiveno rješenje je unutar ℛ1; CPU = ",time[[1]]],Print["Dobiveno rješenje NIJE unutar ℛ1;  CPU = ",time[[1]]]];
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=n j-n+1;
cenDIR[[j]]=Table[zz[[jk+(s-1)]],{s,n}];
,{j,k}];
Print["Initial approximation z0: ",z0=cenDIR];
time1=Timing[z1=Orthogonalize[z0]];
Y0=Table[{Y[[i]],i},{i,m}]; 
(* WKMeans *)
time2=Timing[sol=WKMeansDop[w,Y0,z1]];
sol0=sol[[1]];
z2=sol[[2]];
Print[Style["Solution:",16,Blue]];
Print[Style["F = ",16,Blue],Style[sol[[3]],16,Blue],Style["; Sigma = ",16,Blue],Style[Sigma,16,Blue],"\nThe highest eigenvalues: ",eig[[1]],"\nc = ",Round[sol[[2]],-.001], "; Iteracija: ",sol[[4]],"; CPU = ",time[[1]]];
Print["Ini. Appr.: ",Round[z0,-.0001]];
Print["Orthogonalize: ",Round[z1,-.0001]];
Print["k-means: ",Round[z2,-.0001]];
Print["Time {Ini, Orthogonalize, k-means}: ",{time[[1]],time1[[1]],time2[[1]]}," = ",time[[1]]+time1[[1]]+time2[[1]]];
						(* Konstrukcija particije rješenja *)
PI=Table[{},{j,k}];
Do[
ss=sol0[[j,All,3]];
PI[[j]]=Table[S[[s]],{s,ss}]
,{j,k}];
];
{Y,PI, sol[[2]],sol[[3]],sol[[4]]},
Print[Style["Ifinite expression (1/0) encountered.",14,Red],Style["\nChange the parameter Sigma",14,Blue]]]
]

Graphical presentation of the results 
only for S ⊂ R^2 and S ⊂ R^3

SCRes[sc_]:=Module[{k=Length[sc[[2]]]},
slPI=ListPlot[Table[sc[[2,j]],{j,k}],PlotStyle->{PointSize[.02]},Axes->False,AspectRatio->Automatic,PlotLabel->"Reconstruction"];
Clear[x,y,z];AG=Table[{-1,1},k];
						(* k = 2 *)
If[k==2, OO={0,0};
rj1=Show[ListPlot[Y,PlotStyle->{PointSize->.03} ],
Graphics[{Thick,Gray,Circle[OO,1]}],Graphics[{Gray,Thick,Arrow[{OO,z2[[1]]}],Arrow[{OO,z2[[2]]}]}],AxesOrigin->OO,Ticks->{{-1,1},{-1,1}},AspectRatio->Automatic,
PlotRange->AG,AxesLabel->{x,y}];
rj20=Show[ListPlot[Y,PlotStyle->{PointSize[.03]}  ],
Graphics[{Thick,Gray,Circle[OO,1]}],Graphics[{Gray,Thick,Arrow[{OO,z0[[1]]}],Arrow[{OO,z0[[2]]}]}],AxesOrigin->OO,Ticks->{{-1,1},{-1,1}},AspectRatio->Automatic,
PlotRange->AG,AxesLabel->{x,y},PlotLabel->"Ini. Appr."];
rj21=Show[ListPlot[Y,PlotStyle->{PointSize[.03]}  ],
Graphics[{Thick,Gray,Circle[OO,1]}],Graphics[{Gray,Thick,Arrow[{OO,z1[[1]]}],Arrow[{OO,z1[[2]]}]}],AxesOrigin->OO,Ticks->{{-1,1},{-1,1}},AspectRatio->Automatic,
PlotRange->AG,AxesLabel->{x,y},PlotLabel->"Normalize[]"];
rj22=Show[ListPlot[Y,PlotStyle->{PointSize[.04]}  ],
Graphics[{Thick,Gray,Circle[{0,0},1]}],Graphics[{Gray,Thick,Arrow[{OO,z2[[1]]}],Arrow[{OO,z2[[2]]}]}],AxesOrigin->{0,0},Ticks->{{-1,1},{-1,1}},AspectRatio->Automatic,
PlotRange->AG,AxesLabel->{x,y},PlotLabel->"k-means"];
Print[GraphicsGrid[{{rj20,rj21,rj22}},ImageSize->Large]];
];
							(* k = 3 *)
If[k==3,OO={0,0,0};
rj1=Show[ListPointPlot3D[Y,PlotStyle->{PointSize[.025]} ],
Graphics3D[{Opacity[.2],Sphere[OO,1]}],Graphics3D[{Gray,Thick,Arrow[{OO,z2[[1]]}],Arrow[{OO,z2[[2]]}],Arrow[{OO,z2[[3]]}]}],PlotLabel->"Unit 3-sphere",
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},Ticks->{{1},{1},{1}}];
rj20=Show[ListPointPlot3D[Y,PlotStyle->{PointSize[.025]}],ListPointPlot3D[Y,PlotStyle->{PointSize[.025]}],
Graphics3D[{Opacity[.2],Sphere[OO,1]}],Graphics3D[{Gray,Thick,Arrow[{OO,z0[[1]]}],Arrow[{OO,z0[[2]]}],Arrow[{OO,z0[[3]]}]}],PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},Ticks->{{1},{1},{1}},ImageSize->Medium,PlotLabel->"Ini. Appr."];
rj21=Show[ListPointPlot3D[Y,PlotStyle->{PointSize[.025]}],ListPointPlot3D[z1,PlotStyle->{PointSize[.025]}],
Graphics3D[{Opacity[.2],Sphere[OO,1]}],Graphics3D[{Gray,Thick,Arrow[{OO,z1[[1]]}],Arrow[{OO,z1[[2]]}],Arrow[{OO,z1[[3]]}]}],PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},Ticks->{{1},{1},{1}},ImageSize->Medium,PlotLabel->"Orthogonalize[]"];
rj22=Show[ListPointPlot3D[Y,PlotStyle->{PointSize[.025]}],ListPointPlot3D[z2,PlotStyle->{PointSize[.025]}],
Graphics3D[{Opacity[.2],Sphere[OO,1]}],Graphics3D[{Gray,Thick,Arrow[{OO,z2[[1]]}],Arrow[{OO,z2[[2]]}],Arrow[{OO,z2[[3]]}]}],PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},Ticks->{{1},{1},{1}},ImageSize->Medium,PlotLabel->"k-means"];
Print[GraphicsGrid[{{rj20,rj21,rj22}},ImageSize->Large]]
];
If[k≤3,Print[GraphicsGrid[{{sl0,slPI,rj1}},ImageSize->Large]]]; 
If[k>3,Print[GraphicsGrid[{{sl0,slPI}},ImageSize->Medium]]];
]

Examples

Example-1 (k=2) - IDE

(*SeedRandom[13]*)
Cen={2.5,2.5};sigma=.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
 S=Union[data1,data2]; m=Length[S];k=2; 
Print[Style["|2 kruga| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,AspectRatio->Automatic,ImageSize->Small]
]

Example-2 (k=2) - NE IDE

(*SeedRandom[13]*)
Cen={2.5,2.5};sigma=.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
Cen={6,3};
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
 S=Union[data1,data2]; m=Length[S];k=2;
Print[Style["|2 kruga| = ",16,Blue],Style[m,16,Blue]]; 
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example-3 (k=2) -  IDE!

(*SeedRandom[13]*)
Cen={2.5,2.5};sigma=.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
Cen={10,3};
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
 S=Union[data1,data2]; m=Length[S];k=2; 
Print[Style["|2 kruga| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example-4 (k=2)  - NE IDE

(*SeedRandom[13]*)
cA={0,0}; cB={6,6};broj=10; sigma=0.003;
data1=DataLine[cA,cB,broj,sigma];
cA={6,0}; cB={0,4};
data2=DataLine[cA,cB,broj,sigma];
S=Union[data1,data2]; m=Length[S];k=2; 
Print[Style["|2 pravca| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example-5 (k=2)  -

SeedRandom[13]
Sigma=0.05; (* za .05 opt, za .1 2pts OK, za .15 puno pts OK, za .2- punopts loše *) IT=2;
 Cen={0,4Sqrt[3]/6};sigma=0.002;
data1=DataCircle[Cen,.75,80,sigma,0,2Pi];
cA={-2,0}; cB={2,0};broj=20;
data21=DataLine[cA,cB,broj,sigma];
cA={-2,0}; cB={0,4Sqrt[3]/2};broj=20;
data22=DataLine[cA,cB,broj,sigma];
cA={2,0}; cB={0,4Sqrt[3]/2};broj=20;
data23=DataLine[cA,cB,broj,sigma];
data2=Union[data21,data22,data23];
S=Union[data1,data2]; m=Length[S];k=2; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,PlotRange->All,AspectRatio->Automatic,Axes->None,ImageSize->Small]
]

Example-6 (k=2)  -

SeedRandom[13]
Sigma=0.05; (* za .05 opt, za .1 2pts OK, za .15 puno pts OK, za .2- punopts loše *) IT=2;
 Cen={0,1};sigma=0.002;
data1=DataCircle[Cen,2.25,75,sigma,Pi/2.-.7,Pi/2.+.7];
cA={-1,0}; cB={1,0};broj=20;
data21=DataLine[cA,cB,broj,sigma];
cA={-1,2.5}; cB={1,2.5};broj=20;
data22=DataLine[cA,cB,broj,sigma];
cA={-1,0}; cB={-1,2.5};broj=20;
data23=DataLine[cA,cB,broj,sigma];
cA={1,0}; cB={1,2.5};broj=20;
data24=DataLine[cA,cB,broj,sigma];
data2=Union[data21,data22,data23,data24];
S=Union[data1,data2]; m=Length[S];k=2; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,PlotRange->All,AspectRatio->Automatic,Axes->None,ImageSize->Small]
]

Example-7 (k=3) (broj podataka ovisan o duljini kružnice) - IDE

(*SeedRandom[13]*)
Cen={2.5,2.5};sigma=.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
data3=DataCircle[Cen,1,100,sigma,0,2Pi];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|3 kruga| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,AspectRatio->Automatic,ImageSize->Small]
]

Example-8 (k=3)  (broj podataka ovisan o duljini kružnice) - IDE

SeedRandom[13]
Sigma=0.1;(* .5 ne ide, .15-.35 dobar, .4 dobar, ali puno točaka, .45- puno loš*) IT=5;
Cen={2.5,2.5};sigma=.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
sigma=.1{{1,0},{0,1}};
data3=RandomReal[MultinormalDistribution[Cen,sigma],75];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|3 kruga| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,AspectRatio->Automatic,ImageSize->Small]
]

Example-9 (k=3) - IDE!

(*SeedRandom[13]*)
Cen={4,2.5};sigma=.003;
data1=DataCircle[Cen,3,250,sigma,0,2Pi];
Cen={11,3};
data2=DataCircle[Cen,2.5,150,sigma,0,2Pi];
Cen={-.5,4};
data3=DataCircle[Cen,1,50,sigma,0,2Pi];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|3 kruga| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example-10 (k=3) (jednak broj podataka po kružnicama) - IDE

(*SeedRandom[13]*)
Cen={2.5,2.5};mm1=200;sigma=.003;
data1=DataCircle[Cen,4,mm1,sigma,0,2Pi];
data2=DataCircle[Cen,2.5,mm1,sigma,0,2Pi];
data3=DataCircle[Cen,1,mm1,sigma,0,2Pi];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,AspectRatio->Automatic,ImageSize->Small]
]

Example-11 (k=3)  - NE IDE

(*SeedRandom[13]*)
 Cen={2.5,2.5};sigma=0.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
cA={-1,7}; cB={7.5,-1.5};broj=10;
data3=DataLine[cA,cB,broj,sigma];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,AspectRatio->Automatic,ImageSize->Small]
]

Example-12 (k=3)  - NE IDE

(*SeedRandom[13]*)
cA={0,0}; cB={3,6};broj=10; sigma=0.003;
data1=DataLine[cA,cB,broj,sigma];
cA={6,0}; cB={3,6};
data2=DataLine[cA,cB,broj,sigma];
cA={1,2}; cB={5,2};broj=7;
data3=DataLine[cA,cB,broj,sigma];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|A| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example-13 (k=3)  -  IDE!

(*SeedRandom[13]*)
cA={0,0}; cB={3,5};broj=10; sigma=0.003;
data1=DataLine[cA,cB,broj,sigma];
cA={6,0}; cB={3,4};
data2=DataLine[cA,cB,broj,sigma];
Cen={3,1};  
data3=DataCircle[Cen,1,100,sigma,0,2Pi];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example-14 (k=3)  -  IDE!

(*SeedRandom[13]*)
Sigma=.05; (*ili*)  Sigma=0.1;
Cen={1,2};   sigma=0.003;
data1=DataCircle[Cen,1,100,sigma,0,2Pi];
cA={2.5,1}; cB={2.5,4}; broj=15;
data21=DataLine[cA,cB,broj,sigma];
cA={2.5,2}; cB={3.5,3}; broj=10;
data22=DataLine[cA,cB,broj,sigma];
cA={3,2.5}; cB={4,1}; broj=8;
data23=DataLine[cA,cB,broj,sigma];
data2=Union[data21,data22,data23];
 Cen={5.5,2};
data3=DataCircle[Cen,1,100,sigma,0,2Pi];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|oko| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,PlotRange->All,GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
]

Example -15 (k=6)  -

(*SeedRandom[13]*)
Sigma=0.1;
 (* G *)
Cen={1.5,2};   sigma=0.001;    
data11=DataCircle[Cen,.5,50,sigma,0,Pi];
Cen={1.5,1.5};   
data12=DataCircle[Cen,.5,50,sigma,Pi,2Pi];
cA={1,1.5}; cB={1,2}; broj=20;
data13=DataLine[cA,cB,broj,sigma];
cA={1.5,1.5}; cB={2,1.5}; broj=20;
data14=DataLine[cA,cB,broj,sigma];
data1=Union[data11,data12,data13,data14];
(* A *)
Cen={3,2};          sigma=0.001;                                               
data21=DataCircle[Cen,.5,50,sigma,0,Pi];
cA={2.5,1}; cB={2.5,2}; broj=30;
data22=DataLine[cA,cB,broj,sigma];
cA={3.5,1}; cB={3.5,2}; broj=30;
data23=DataLine[cA,cB,broj,sigma];
cA={2.5,1.5}; cB={3.5,1.5}; broj=20;
data24=DataLine[cA,cB,broj,sigma];
data2=Union[data21,data22,data23,data24];
(* U *)
Cen={4.5,1.5};   sigma=0.001;     
data31=DataCircle[Cen,.5,50,sigma,Pi,2Pi];
cA={4,1.5}; cB={4,2.5}; broj=30;
data32=DataLine[cA,cB,broj,sigma];
cA={5,1.5}; cB={5,2.5}; broj=30;
data33=DataLine[cA,cB,broj,sigma];
data3=Union[data31,data32,data33];
(* S *)
Cen={6,1.375};   sigma=0.001;     
data41=DataCircle[Cen,.375,50,sigma,-Pi/2,Pi/2];
Cen={6,2.125};  
data42=DataCircle[Cen,.375,50,sigma,Pi/2,3Pi/2];
cA={6,2.5}; cB={6.5,2.5}; broj=20;
data43=DataLine[cA,cB,broj,sigma];
cA={6,1}; cB={5.5,1}; broj=20;
data44=DataLine[cA,cB,broj,sigma];
data4=Union[data41,data42,data43,data44];
(* S *)
Cen={7.5,1.375};   sigma=0.001;     
data51=DataCircle[Cen,.375,50,sigma,-Pi/2,Pi/2];
Cen={7.5,2.125};  
data52=DataCircle[Cen,.375,50,sigma,Pi/2,3Pi/2];
cA={7.5,2.5}; cB={8,2.5}; broj=20;
data53=DataLine[cA,cB,broj,sigma];
cA={7.5,1}; cB={7,1}; broj=20;
data54=DataLine[cA,cB,broj,sigma];
data5=Union[data51,data52,data53,data54];
(*  LINE *)
cA={0,2}; cB={1,.5}; broj=10; sigma=0.001; 
data61=DataLine[cA,cB,broj,sigma];
cA={1,.5}; cB={8,.5}; broj=20;
data62=DataLine[cA,cB,broj,sigma];
data6=Union[data61,data62];

S=Union[data1,data2,data3,data4,data5,data6]; m=Length[S];k=6; 
Print[Style["|Gauss| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.01]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.01]}];
sl3=ListPlot[data3,PlotStyle->{Red,PointSize[.01]}];
sl4=ListPlot[data4,PlotStyle->{Cyan,PointSize[.01]}];
sl5=ListPlot[data5,PlotStyle->{Magenta,PointSize[.01]}];
sl6=ListPlot[data6,PlotStyle->{Gray,PointSize[.01]}];
sl0=Show[sl1,sl2,sl3,sl4,sl5,sl6,PlotRange->All,AspectRatio->Automatic,Axes->False,ImageSize->Medium]
]

Example-16 (k=4)  -  IDE!

SeedRandom[13]
Sigma=0.1; IT=10;   (* optimal choice of Initial approximation - NEKI PUTA NE IDE *)
a=2;b=1; sigma=0.03;
data1=DataEll[a,b,200,sigma];
Cen={1,0}; sigma=.01{{2,1},{1,2}};
data2=RandomReal[MultinormalDistribution[Cen,sigma],75];
cA={-2,1}; cB={0,2}; broj=20; sigma=0.001;
data31=DataLine[cA,cB,broj,sigma];
cA={2,1}; cB={0,2}; broj=20; sigma=0.001;
data32=DataLine[cA,cB,broj,sigma];
data3=Union[data31,data32];
Cen={0,3};   sigma=0.001;     
data4=DataCircle[Cen,4.5,100,sigma,21.5Pi/16,26.5Pi/16];

S=Union[data1,data2,data3,data4]; m=Length[S];k=4; 
Print[Style["|Ellipse| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]},Axes->None];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Cyan,PointSize[.02]}];
sl4=ListPlot[data4,PlotStyle->{Magenta,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,sl4,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example - 17 (k = 5)  -  IDE!

SeedRandom[13]
Sigma=0.1; IT=10;   (*   *)
Cen={0,4.5};   sigma=0.001;     
data1=DataCircle[Cen,2,100,sigma,Pi+.2,2Pi-.2];
Cen={-2.5,3.5};   sigma=0.001;     
data2=DataCircle[Cen,2,100,sigma,0,1.3Pi];
Cen={2.5,3.5};   sigma=0.001;    
data3=DataCircle[Cen,2,100,sigma,-.7Pi/2,Pi];
Cen={3,.5};   sigma=0.001;    
data4=DataCircle[Cen,2,100,sigma,Pi/2,Pi];
Cen={-3,.5};   sigma=0.001;    
data5=DataCircle[Cen,2,100,sigma,0,Pi/2];


S=Union[data1,data2,data3,data4,data5]; m=Length[S];k=5; 
Print[Style["|Vijenac| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]},Axes->None];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Red,PointSize[.02]}];
sl4=ListPlot[data4,PlotStyle->{Cyan,PointSize[.02]}];
sl5=ListPlot[data5,PlotStyle->{Magenta,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,sl4,sl5,PlotRange->All, AspectRatio->Automatic,ImageSize->Small]
]

Examples chosen

Example-5 (k=2)  -

SeedRandom[13]
Sigma=0.05; (* za .05 opt, za .1 2pts OK, za .15 puno pts OK, za .2- punopts loše *) IT=2;
(* ,{sig,0.05,0.45,0.025} *)
 Cen={0,4Sqrt[3]/6};sigma=0.002;
data1=DataCircle[Cen,.75,80,sigma,0,2Pi];
cA={-2,0}; cB={2,0};broj=20;
data21=DataLine[cA,cB,broj,sigma];
cA={-2,0}; cB={0,4Sqrt[3]/2};broj=20;
data22=DataLine[cA,cB,broj,sigma];
cA={2,0}; cB={0,4Sqrt[3]/2};broj=20;
data23=DataLine[cA,cB,broj,sigma];
data2=Union[data21,data22,data23];
S=Union[data1,data2]; m=Length[S];k=2; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl0=Show[sl1,sl2,PlotRange->All,AspectRatio->Automatic,Axes->None,PlotLabel->"Dataset",ImageSize->Small]
]

Example-8 (k=3)  (broj podataka ovisan o duljini kružnice) - IDE

SeedRandom[13] 
Sigma=0.1; IT=5;(* .5 ne ide, .15-.35 dobar, .4 dobar, ali puno točaka, .45- puno loš;
{sig,0.05,0.45,0.025} *)
Cen={2.5,2.5};sigma=.003;
data1=DataCircle[Cen,4,300,sigma,0,2Pi];
data2=DataCircle[Cen,2.5,200,sigma,0,2Pi];
sigma=.1{{1,0},{0,1}};
data3=RandomReal[MultinormalDistribution[Cen,sigma],75];
S=Union[data1,data2,data3]; m=Length[S];k=3; 
Print[Style["|S| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Gray,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,AspectRatio->Automatic,Axes->None,PlotLabel->"Data set",ImageSize->Small]
]

Example -15 (k=6)  -

(*SeedRandom[13]*)
Sigma=0.1;
 (* G *)
Cen={1.5,2};   sigma=0.001;    
data11=DataCircle[Cen,.5,50,sigma,0,Pi];
Cen={1.5,1.5};   
data12=DataCircle[Cen,.5,50,sigma,Pi,2Pi];
cA={1,1.5}; cB={1,2}; broj=20;
data13=DataLine[cA,cB,broj,sigma];
cA={1.5,1.5}; cB={2,1.5}; broj=20;
data14=DataLine[cA,cB,broj,sigma];
data1=Union[data11,data12,data13,data14];
(* A *)
Cen={3,2};          sigma=0.001;                                               
data21=DataCircle[Cen,.5,50,sigma,0,Pi];
cA={2.5,1}; cB={2.5,2}; broj=30;
data22=DataLine[cA,cB,broj,sigma];
cA={3.5,1}; cB={3.5,2}; broj=30;
data23=DataLine[cA,cB,broj,sigma];
cA={2.5,1.5}; cB={3.5,1.5}; broj=20;
data24=DataLine[cA,cB,broj,sigma];
data2=Union[data21,data22,data23,data24];
(* U *)
Cen={4.5,1.5};   sigma=0.001;     
data31=DataCircle[Cen,.5,50,sigma,Pi,2Pi];
cA={4,1.5}; cB={4,2.5}; broj=30;
data32=DataLine[cA,cB,broj,sigma];
cA={5,1.5}; cB={5,2.5}; broj=30;
data33=DataLine[cA,cB,broj,sigma];
data3=Union[data31,data32,data33];
(* S *)
Cen={6,1.375};   sigma=0.001;     
data41=DataCircle[Cen,.375,50,sigma,-Pi/2,Pi/2];
Cen={6,2.125};  
data42=DataCircle[Cen,.375,50,sigma,Pi/2,3Pi/2];
cA={6,2.5}; cB={6.5,2.5}; broj=20;
data43=DataLine[cA,cB,broj,sigma];
cA={6,1}; cB={5.5,1}; broj=20;
data44=DataLine[cA,cB,broj,sigma];
data4=Union[data41,data42,data43,data44];
(* S *)
Cen={7.5,1.375};   sigma=0.001;     
data51=DataCircle[Cen,.375,50,sigma,-Pi/2,Pi/2];
Cen={7.5,2.125};  
data52=DataCircle[Cen,.375,50,sigma,Pi/2,3Pi/2];
cA={7.5,2.5}; cB={8,2.5}; broj=20;
data53=DataLine[cA,cB,broj,sigma];
cA={7.5,1}; cB={7,1}; broj=20;
data54=DataLine[cA,cB,broj,sigma];
data5=Union[data51,data52,data53,data54];
(*  LINE *)
cA={0,2}; cB={1,.5}; broj=10; sigma=0.001; 
data61=DataLine[cA,cB,broj,sigma];
cA={1,.5}; cB={8,.5}; broj=20;
data62=DataLine[cA,cB,broj,sigma];
data6=Union[data61,data62];

S=Union[data1,data2,data3,data4,data5,data6]; m=Length[S];k=6; 
Print[Style["|Gauss| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.01]}];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.01]}];
sl3=ListPlot[data3,PlotStyle->{Red,PointSize[.01]}];
sl4=ListPlot[data4,PlotStyle->{Cyan,PointSize[.01]}];
sl5=ListPlot[data5,PlotStyle->{Magenta,PointSize[.01]}];
sl6=ListPlot[data6,PlotStyle->{Gray,PointSize[.01]}];
sl0=Show[sl1,sl2,sl3,sl4,sl5,sl6,PlotRange->All,AspectRatio->Automatic,Axes->False,ImageSize->Medium]
]

Example-16 (k=4)  -  IDE!

SeedRandom[13]
Sigma=0.1; IT=10;   (* optimal choice of Initial approximation - NEKI PUTA NE IDE *)
a=2;b=1; sigma=0.03;
data1=DataEll[a,b,200,sigma];
Cen={1,0}; sigma=.01{{2,1},{1,2}};
data2=RandomReal[MultinormalDistribution[Cen,sigma],75];
cA={-2,1}; cB={0,2}; broj=20; sigma=0.001;
data31=DataLine[cA,cB,broj,sigma];
cA={2,1}; cB={0,2}; broj=20; sigma=0.001;
data32=DataLine[cA,cB,broj,sigma];
data3=Union[data31,data32];
Cen={0,3};   sigma=0.001;     
data4=DataCircle[Cen,4.5,100,sigma,21.5Pi/16,26.5Pi/16];

S=Union[data1,data2,data3,data4]; m=Length[S];k=4; 
Print[Style["|Ellipse| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]},Axes->None];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Cyan,PointSize[.02]}];
sl4=ListPlot[data4,PlotStyle->{Magenta,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,sl4,PlotRange->All,AspectRatio->Automatic,ImageSize->Small]
]

Example - 17 (k = 5)  -  IDE!

SeedRandom[13]
Sigma=0.1; IT=10;   (*   *)
Cen={0,4.5};   sigma=0.001;     
data1=DataCircle[Cen,2,100,sigma,Pi+.2,2Pi-.2];
Cen={-2.5,3.5};   sigma=0.001;     
data2=DataCircle[Cen,2,100,sigma,0,1.3Pi];
Cen={2.5,3.5};   sigma=0.001;    
data3=DataCircle[Cen,2,100,sigma,-.7Pi/2,Pi];
Cen={3,.5};   sigma=0.001;    
data4=DataCircle[Cen,2,100,sigma,Pi/2,Pi];
Cen={-3,.5};   sigma=0.001;    
data5=DataCircle[Cen,2,100,sigma,0,Pi/2];

S=Union[data1,data2,data3,data4,data5]; m=Length[S];k=5; 
Print[Style["|Vijenac| = ",16,Blue],Style[m,16,Blue]];
If[Length[data1[[1]]]==2,
sl1=ListPlot[data1,PlotStyle->{Blue,PointSize[.02]},Axes->None];
sl2=ListPlot[data2,PlotStyle->{Orange,PointSize[.02]}];
sl3=ListPlot[data3,PlotStyle->{Red,PointSize[.02]}];
sl4=ListPlot[data4,PlotStyle->{Cyan,PointSize[.02]}];
sl5=ListPlot[data5,PlotStyle->{Magenta,PointSize[.02]}];
sl0=Show[sl1,sl2,sl3,sl4,sl5,PlotRange->All, AspectRatio->Automatic,ImageSize->Small]
]

Processing

S: (* a linearly non -separable dataset *);
k: number of clusters;
Sigma: parameter [0.1, 5];
Alg: Spectar[1], SpectarLCDIR[2];
IT: number of iterations for SpectarLCDIR

(* SCAlg[S_,k_,Sigma_,Alg_,IT_] *)
sc=SCAlg[S,6,0.1,1,5];
SCRes[sc]