Moduli

Needs["ComputationalGeometry`"]

Transformacije

T[x_] :=Module[{n=Length[x],dd,MD,vu},
dd = Table[AG[[i, 2]] - AG[[i, 1]], {i, n}];
MD = DiagonalMatrix[1/dd]; 
vu = Table[AG[[i,1]], {i, n}];
MD.(x - vu)
]
T1[x_] :=Module[{n=Length[x],dd,MD1,vu},
dd = Table[AG[[i, 2]] - AG[[i, 1]], {i, n}];
MD1 = DiagonalMatrix[dd];
vu = Table[AG[[i,1]], {i, n}];
MD1.x + vu
]

Hausdorff distance between two sets

Input: finite sets A,B ⊂ (ℝ^n); 
Output: Hausdorff distance;
Submodules: none needed

Hdist[a_,b_]:=Module[{m,n,a1,b1},
If[Length[a]>Length[b],a1=b;b1=a,a1=a;b1=b];
m=Length[a1];n=Length[b1];
Max[Max[Table[Min[Table[Norm[a1[[i]]-b1[[j]]],{j,n}]],{i,m}]],Max[Table[Min[Table[Norm[a1[[i]]-b1[[j]]],{i,m}]],{j,n}]]]]

Weighted k -  Means

WKMeans[A_,w_,z_]:=Module[{m=Length[A],k=Length[z],n=Length[A[[1]]],c0,Pod,imin,PI,IT=5,It=0},
Pod=Table[{w[[i]],A[[i]]},{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]; 
						(* 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"]];
Chop[Norm[c-c0]]≠0&& It<IT,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];  
               ];
{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[A_,w_,z_]:=Module[{},
Sum[w[[i]]Min[Table[Norm[z[[j]]-A[[i]]]^2,{j,Length[z]}]],{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]} ]]
];

Rand and Jaccard index and the confusion matrix

Input: two partitions (Π^((1))), (Π^((2))) of the set 𝒜 ⊂ (ℝ^n); 
Output: Rand index R, Jaccard index J, confusion matrix S, CPU-time;
Submodules: none needed

Rand[PI1_,PI2_]:=Module[{m=Length[Flatten[PI1]],S,S1=Table[0,{i,Length[PI1]}],S2=Table[0,{j,Length[PI2]}],a,b,c,d},
time=Timing[
(* Confusion matrix *)
S=Table[0,{r,Length[PI1]},{r,Length[PI2]}];
Do[
Do[
S[[r,s]]=Length[Intersection[PI1[[r]],PI2[[s]]]]
,{r,Length[PI1]}]
,{s,Length[PI2]}];
(* Numbers a, b, c, d *)
TT=Sum[S[[i,j]]^2,{i,Length[PI1]},{j,Length[PI2]}];
Do[
S1[[i]]=Sum[S[[i,j]],{j,Length[PI2]}]
,{i,Length[PI1]}];
Do[
S2[[j]]=Sum[S[[i,j]],{i,Length[PI1]}]
,{j,Length[PI2]}];
a=.5(TT-m);
b=.5(Sum[S1[[i]]^2,{i,Length[PI1]}]-TT);
c=.5(Sum[S2[[j]]^2,{j,Length[PI2]}]-TT);
d=.5(TT+m^2-Sum[S1[[i]]^2,{i,Length[PI1]}]-Sum[S2[[j]]^2,{j,Length[PI2]}]);
];
{(a+d)/Binomial[m,2],a/(a+b+c),S,time[[1]]}
]

DIRECT

Modul for solving the global optimization problem for a Lipschitz continuous function. 
Mathematica codes written by I.Vazler;
Input: ordinary or compiled function, Domain , {Eps->10^(-4),MaxIter->20,DMin->10.^(-2)},
Output: optimal variables values, optimal objective function value , CPU-time, ;
Submodules: none needed

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,dtest,anim={}},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
DIRECT::maxiter="The maximal number `1` of iterations has been reached.";
n=Length[dom];Nf=0;
iter=0;
IT=0;
(* Scaling the function onto the unit hypercube *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* Initial centers and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]}; Nf=Nf+1;
TMin={c[[1]],fc[[1]]};
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: Searching for the potentially optimal hyperrectangle *)
dmaxrazni=Union[dmax];
(* all data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* data with the lawest f-value grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* checking the rest of the requirements *)
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 the potentially optimal directions *)
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]}];
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}];
(* Printing ordinal number of the iteration and number of potentially optimal rectangles 
Print["Iter: ",iter," " ,Length[pot2]]; *)
(* Printing iterations
Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* figure for the 2D case *)
(*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]}];
IT=iter;
(* back to the original domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* printing the figure in the 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
IT=iter;
Return[TMin];
]

SymDIRECT

Modul for solving the global optimization problem for a simmetric Lipschitz continuous function. 
Mathematica codes written by I.Vazler;
Input: ordinary or compiled function, Domain ,
(*{Eps->10^(-4),MaxIter->20,DMin->10.^(-2)} *);
Output: optimal variables values, optimal objective function value , CPU-time, ;
Submodules: none needed

(* For a given center and length of the side of the hyperrectangle the function test returns:
  1 If the whole hyperrectangle is located in the region \Delta
0 If a hyperrectangle lies in the region \Delta only partially
-1 If a hyperrectangle appears outside the region \Delta (or if a hyperrectangle has only one point in \Delta)*)
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,dtest,anim={},unutra,brisi},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
SymDIRECT::maxiter="The maximal number `1` of iterations has been reached.";
iter=0; Nf=0;
n=Length[dom];
(* 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={test[c[[1]],d[[1]]]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};   Nf=Nf+1;
TMin={c[[1]],fc[[1]]};
(* 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 *)
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];
,
(* 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 *)
temp=Position[novifc,Min[novifc]][[1,1]];
If[novifc[[temp]]<TMin[[2]],TMin=N[{novic[[temp]],novifc[[temp]]},10]];
,{pot,pot2}];
(* 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];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[SymDIRECT::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];
]

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="The maximal number `1` of iterations has been reached.";
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 - LP

Modul koji provjerava nalazi li se  rectangle R[c,dd];
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

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 ;
    ];
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;
Ovo smo isključili (iako se primjerice za n k=13 CPU-time skrati se 10 puta; 
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
    ];
tt
]

Example : 4 - 6 - 8 - 10 clusters in R^2

4 clusters

SeedRandom[13];mc=100;
c={{.4,.3},{.25,.55},{.75,.3},{.6,.75}}; k=Length[c]; n=Length[c[[1]]];
AG={{0,1},{0,1}};
kov=Table[IdentityMatrix[2].005,{i,k}];
podaci=Table[RandomReal[MultinormalDistribution[c[[i]],kov[[i]]],mc],{i,k}];
B=Flatten[podaci,1];
m=Length[B];
w=Table[1,{i,m}];
sol=WKMeans[B,w,c];
c=sol[[2]];
Print["Data set A: c0 = ",sol[[2]],"; F0 = ",F0=sol[[3]]];
PI0=Table[sol[[1,j,All,2]],{j,k}];
Print["Nr. of elements in clusters: ",Table[Length[PI0[[r]]],{r,k}]];
slc=ListPlot[c,PlotMarkers->{"✶",18},PlotStyle->{Brown,PointSize[.04],Opacity[.95]}];
slPI0=ListPlot[PI0,AspectRatio->Automatic,PlotStyle->{PointSize[.02]},PlotRange->AG,ImageSize->Small];
sl1=Show[slPI0,slc,GridLines->Automatic,ImageSize->Small];
slB=ListPlot[B,AspectRatio->Automatic,PlotStyle->{PointSize[.02]}];
sl2=Show[slB,slc,PlotRange->{{0,1},{0,1}},GridLines->Automatic,ImageSize->Small];
Print[Show[GraphicsGrid[{{sl1,sl2}}],ImageSize->Medium]]
					(* Ograničenja *)
kn= k n;  
Nf=LP=LP1=0;
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}];
Print["[A1;b1] = ",{MatrixForm[A1],MatrixForm[b1]}]
Print["{m1,n1} = ",{m1,n1}=Dimensions[A1],"; pts = ",pts]

Execution

ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4,x5,x6,x7,x8}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4,x5,x6,x7,x8}]
dom=Table[{0,1},{kn}]; Nf=0;
							(* COMPILE *)
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];
tDIR=Timing[min=LinConstrDIRECTe[FFCompile,dom,DMin->10.^(-3),MaxIter->10, Eps->10^(-4)]]
If[RegionMember[ℛ1,min[[1]]]==True,Print["Obtained solution is INSIDE of ℛ1"],Print["Obtained solution is NOT INSIDE of ℛ1"]];
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Calculated centers in [0,1]^kn: ",cenDIR,";  Hdist = ",
Hdir=Hdist[c,cenDIR],"; Nf = ",Nf,"; F = ",min[[2]]];
Print[Style["Strategy 1 is used",Blue,18]];
						(* WKmeans *)
twk=Timing[sol=WKMeans[B,w,cenDIR]];
slcN=ListPlot[sol[[2]],PlotStyle->{Red,PointSize[.02],Opacity[.9]}];
PI1=Table[sol[[1,j,All,2]],{j,k}];
Show[slB,slc,slcN,GridLines->Automatic,ImageSize->Small]
Print["WKMeans centers in [0,1]^kn: ",sol[[2]],"\nHdist = ",
Hwk=Hdist[c,sol[[2]]],"; F = ",sol[[3]],"; It = ",sol[[4]],"; CPUtotal=",tDIR[[1]]+twk[[1]],"; Rand = ",rand=Rand[PI0,PI1],"; CPU = ",time[[1]],"\nConf.Mat. = ",rand[[3]]];
Print["Alldata (k,m,F0,F1,Hdir,cpu,it,F2,Hwk,Rand,cpu): ",{k,m,Round[F0,.01],Round[min[[2]],.01],Round[Hdir,.01],Round[tDIR[[1]],.01],sol[[4]],Round[sol[[3]],.01],Round[Hwk,.01],Round[rand[[1]],.01],Round[twk[[1]],.01]}]

DIRECT - execution

time=Timing[min=DIRECT[FFCompile,dom,DMin->10.^(-4),MaxIter->200, Eps->10^(-4)]]
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Izračunati centri u [0,1]^kn: ",cenDIR,";  Hdist = ",
Hdist[c,cenDIR],"; F = ",min[[2]]];
Print["Alldata (iter,f,Nf,T): ",{iter,min[[2]],Nf,time[[1]]}];

6 clusters

SeedRandom[13];mc=100;
c={{.8,.3},{.65,.5},{.5,.6},{.3,.8},{.4,.2},{.2,.4}};  k=Length[c]; n=Length[c[[1]]];
AG={{0,1},{0,1}};
kov=Table[IdentityMatrix[2].005,{i,k}];
podaci=Table[RandomReal[MultinormalDistribution[c[[i]],kov[[i]]],mc],{i,k}];
B=Flatten[podaci,1];
m=Length[B];
w=Table[1,{i,m}];
sol=WKMeans[B,w,c];
c=sol[[2]];
Print["Data set A: c0 = ",sol[[2]],"; F0 = ",F0=sol[[3]]];
PI0=Table[sol[[1,j,All,2]],{j,k}];
Print["Nr. of elements in clusters: ",Table[Length[PI0[[r]]],{r,k}]];
slc=ListPlot[c,PlotMarkers->{"✶",18},PlotStyle->{Brown,PointSize[.04],Opacity[.95]}];
slPI0=ListPlot[PI0,AspectRatio->Automatic,PlotStyle->{PointSize[.02]},PlotRange->AG,ImageSize->Small];
sl1=Show[slPI0,slc,GridLines->Automatic,ImageSize->Small];
slB=ListPlot[B,AspectRatio->Automatic,PlotStyle->{PointSize[.02]}];
sl2=Show[slB,slc,PlotRange->{{0,1},{0,1}},GridLines->Automatic,ImageSize->Small];
Print[Show[GraphicsGrid[{{sl1,sl2}}],ImageSize->Medium]]
					(* Ograničenja *)
kn= k n;  
Nf=LP=LP1=0;
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}];
Print["[A1;b1] = ",{MatrixForm[A1],MatrixForm[b1]}]
Print["{m1,n1} = ",{m1,n1}=Dimensions[A1],"; pts = ",pts]

Execution

ℛ1=ImplicitRegion[And@@Table[A1[[i]].{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12}≤  
   b1[[i]],{i,Length[b1]}],{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12}]
dom=Table[{0,1},{kn}]; Nf=0;
							(* COMPILE *)
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];
tDIR=Timing[min=LinConstrDIRECTe[FFCompile,dom,DMin->10.^(-3),MaxIter->10, Eps->10^(-4)]]
If[RegionMember[ℛ1,min[[1]]]==True,Print["Obtained solution is INSIDE of ℛ1; |ctest| = ",Length[ctest]],Print["Obtained solution is NOT INSIDE of ℛ1"]];
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Calculated centers in [0,1]^kn: ",cenDIR,";  Hdist = ",
Hdir=Hdist[c,cenDIR],"; Nf = ",Nf,"; F = ",min[[2]]];
Print[Style["Strategy 1 is used",Blue,18]];
						(* WKmeans *)
twk=Timing[sol=WKMeans[B,w,cenDIR]];
slcN=ListPlot[sol[[2]],PlotStyle->{Red,PointSize[.02],Opacity[.9]}];
PI1=Table[sol[[1,j,All,2]],{j,k}];
Show[slB,slc,slcN,GridLines->Automatic,ImageSize->Small]
Print["WKMeans centers in [0,1]^kn: ",sol[[2]],"\nHdist = ",
Hwk=Hdist[c,sol[[2]]],"; F = ",sol[[3]],"; It = ",sol[[4]],"; CPUtotal=",tDIR[[1]]+twk[[1]],"; Rand = ",rand=Rand[PI0,PI1],"; CPU = ",time[[1]],"\nConf.Mat. = ",rand[[3]]];
Print["Alldata (k,kn,m,F0,F1,Hdir,cpu,it,F2,Hwk,Rand,cpu): ",{k,kn,m,Round[F0,.01],Round[min[[2]],.01],Round[Hdir,.01],Round[tDIR[[1]],.01],sol[[4]],Round[sol[[3]],.01],Round[Hwk,.01],Round[rand[[1]],.01],Round[twk[[1]],.01]}]

DIRECT - execution

time=Timing[min=DIRECT[FFCompile,dom,DMin->10.^(-4),MaxIter->200, Eps->10^(-4)]]
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Izračunati centri u [0,1]^kn: ",cenDIR,"\nHdist = ",
Hdist[c1,cenDIR],"; F = ",min[[2]],"; |ctest| = ",Length[ctest]];
Print["Alldata (iter,f,Nf,T): ",{iter,min[[2]],Nf,time[[1]]}];

8 clusters

SeedRandom[13];mc=100;
c={{.7,.2},{.6,.4},{.5,.6},{.25,.75},{.35,.2},{.18,.4},{.83,.6},{.6,.8}};  k=Length[c]; n=Length[c[[1]]];
AG={{0,1},{0,1}};
kov=Table[IdentityMatrix[2].005,{i,k}];
podaci=Table[RandomReal[MultinormalDistribution[c[[i]],kov[[i]]],mc],{i,k}];
B=Flatten[podaci,1];
m=Length[B];
w=Table[1,{i,m}];
sol=WKMeans[B,w,c];
c=sol[[2]];
Print["Data set A: c0 = ",sol[[2]],"; F0 = ",F0=sol[[3]]];
PI0=Table[sol[[1,j,All,2]],{j,k}];
Print["Nr. of elements in clusters: ",Table[Length[PI0[[r]]],{r,k}]];
slc=ListPlot[c,PlotMarkers->{"✶",18},PlotStyle->{Brown,PointSize[.04],Opacity[.95]}];
slPI0=ListPlot[PI0,AspectRatio->Automatic,PlotStyle->{PointSize[.02]},PlotRange->AG,ImageSize->Small];
sl1=Show[slPI0,slc,GridLines->Automatic,ImageSize->Small];
slB=ListPlot[B,AspectRatio->Automatic,PlotStyle->{PointSize[.02]}];
sl2=Show[slB,slc,PlotRange->{{0,1},{0,1}},GridLines->Automatic,ImageSize->Small];
Print[Show[GraphicsGrid[{{sl1,sl2}}],ImageSize->Medium]]
					(* Ograničenja *)
kn= k n;  
Nf=LP=LP1=0;
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}];
Print["[A1;b1] = ",{MatrixForm[A1],MatrixForm[b1]}]
Print["{m1,n1} = ",{m1,n1}=Dimensions[A1],"; pts = ",pts]

Execution

ℛ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}]
dom=Table[{0,1},{kn}]; Nf=0;
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];
tDIR=Timing[min=LinConstrDIRECTe[FFCompile,dom,DMin->10.^(-3),MaxIter->10, Eps->10^(-4)]]
If[RegionMember[ℛ1,min[[1]]]==True,Print["Obtained solution is INSIDE of ℛ1; |ctest| = ",Length[ctest]],Print["Obtained solution is NOT INSIDE of ℛ1"]];
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Calculated centers in [0,1]^kn: ",cenDIR,";  Hdist = ",
Hdir=Hdist[c,cenDIR],"; Nf = ",Nf,"; F = ",min[[2]]];
Print[Style["Strategy 1 is used",Blue,18]];
						(* WKmeans *)
twk=Timing[sol=WKMeans[B,w,cenDIR]];
slcN=ListPlot[sol[[2]],PlotStyle->{Red,PointSize[.02],Opacity[.9]}];
PI1=Table[sol[[1,j,All,2]],{j,k}];
Show[slB,slc,slcN,GridLines->Automatic,ImageSize->Small]
Print["WKMeans centers in [0,1]^kn: ",sol[[2]],"\nHdist = ",
Hwk=Hdist[c,sol[[2]]],"; F = ",sol[[3]],"; It = ",sol[[4]],"; CPUtotal=",tDIR[[1]]+twk[[1]],"; Rand = ",rand=Rand[PI0,PI1],"; CPU = ",time[[1]],"\nConf.Mat. = ",rand[[3]]];
Print["Alldata (k,kn,m,F0,F1,Hdir,cpu,it,F2,Hwk,Rand,cpu): ",{k,kn,m,Round[F0,.01],Round[min[[2]],.01],Round[Hdir,.01],Round[tDIR[[1]],.01],sol[[4]],Round[sol[[3]],.01],Round[Hwk,.01],Round[rand[[1]],.01],Round[twk[[1]],.01]}]

DIRECT - execution

time=Timing[min=DIRECT[FFCompile,dom,DMin->10.^(-4),MaxIter->200, Eps->10^(-4)]]
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Izračunati centri u [0,1]^kn: ",cenDIR,";  Hdist = ",
Hdist[c1,cenDIR],"; F = ",min[[2]]];
Print["Alldata (iter,f,Nf,T): ",{iter,min[[2]],Nf,time[[1]]}];

Za ovih 83 iteracije trebalo mu je 3 dana!!!!!!!!!!!!!

10 clusters

SeedRandom[13];mc=100;
c={{.78,.82},{.83,.55},{.55,.75},{.23,.82},{.5,.5},{.75,.2},{.6,.37},{.3,.6},{.35,.2},{.18,.35}};  k=Length[c]; n=Length[c[[1]]];
AG={{0,1},{0,1}};
kov=Table[IdentityMatrix[2].005,{i,k}];
podaci=Table[RandomReal[MultinormalDistribution[c[[i]],kov[[i]]],mc],{i,k}];
B=Flatten[podaci,1];
m=Length[B];
w=Table[1,{i,m}];
sol=WKMeans[B,w,c];
c=sol[[2]];
Print["Data set A: c0 = ",sol[[2]],"; F0 = ",F0=sol[[3]]];
PI0=Table[sol[[1,j,All,2]],{j,k}];
Print["Nr. of elements in clusters: ",Table[Length[PI0[[r]]],{r,k}]];
slc=ListPlot[c,PlotMarkers->{"✶",15},PlotStyle->{Brown,PointSize[.04],Opacity[.95]}];
slPI0=ListPlot[PI0,AspectRatio->Automatic,PlotStyle->{PointSize[.02]},PlotRange->AG,ImageSize->Small];
(*R={1.75,3.75,5,6.5,8}/4; 
tt=Table[ContourPlot[Norm[{x,y},1]==R[[j]],{x,0,1.1},{y,0,1.1},ContourStyle->Red],{j,Length[R]}];*)
sl1=Show[slPI0,slc,GridLines->Automatic,ImageSize->Small];
slB=ListPlot[B,AspectRatio->Automatic,PlotStyle->{PointSize[.02]}];
sl2=Show[slB,slc,PlotRange->{{0,1},{0,1}},GridLines->Automatic,ImageSize->Small];
Print[Show[GraphicsGrid[{{sl1,sl2}}],ImageSize->Medium]]
					(* Ograničenja *)
kn= k n;  
Nf=LP=LP1=0;
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}];
Print["[A1;b1] = ",{MatrixForm[A1],MatrixForm[b1]}]
Print["{m1,n1} = ",{m1,n1}=Dimensions[A1],"; pts = ",pts]

Execution

ℛ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}≤  
   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}]
dom=Table[{0,1},{kn}]; Nf=0;
			(* COMPILE *)			cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];
tDIR=Timing[min=LinConstrDIRECTe[FFCompile,dom,DMin->10.^(-3),MaxIter->10, Eps->10^(-4)]]
If[RegionMember[ℛ1,min[[1]]]==True,Print["Obtained solution is INSIDE of ℛ1; |ctest| = ",Length[ctest]],Print["Obtained solution is NOT INSIDE of ℛ1"]];
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=2j-1;
cenDIR[[j]]={zz[[jk]],zz[[jk+1]]};
,{j,k}];
Print["Calculated centers in [0,1]^kn: ",cenDIR,";  Hdist = ",
Hdir=Hdist[c,cenDIR],"; Nf = ",Nf,"; F = ",min[[2]]];
Print[Style["Strategy 1 is used",Blue,18]];
						(* WKmeans *)
twk=Timing[sol=WKMeans[B,w,cenDIR]];
slcN=ListPlot[sol[[2]],PlotStyle->{Red,PointSize[.02],Opacity[.9]}];
PI1=Table[sol[[1,j,All,2]],{j,Length[sol[[2]]]}];
Show[slB,slc,slcN,GridLines->Automatic,ImageSize->Small]
Print["WKMeans centers in [0,1]^kn: ",sol[[2]],"\nHdist = ",
Hwk=Hdist[c,sol[[2]]],"; F = ",sol[[3]],"; It = ",sol[[4]],"; CPUtotal=",tDIR[[1]]+twk[[1]],"; Rand = ",rand=Rand[PI0,PI1],"; CPU = ",time[[1]],"\nConf.Mat. = ",rand[[3]]];
Print["Alldata (k,kn,m,F0,F1,Hdir,cpu,it,F2,Hwk,Rand,cpu): ",{k,kn,m,Round[F0,.01],Round[min[[2]],.01],Round[Hdir,.01],Round[tDIR[[1]],.01],sol[[4]],Round[sol[[3]],.01],Round[Hwk,.01],Round[rand[[1]],.01],Round[twk[[1]],.01]}]