Oval Recognizing

(* Oval Recognizing is a complete Mathematica-program for oval recognizing
by using a (Subscript[ℓ, 1]) or LS-distance like function. All modules needed are listed in the preamble. Najprije treba učitati bazu podataka "E-ColiA.txt" i redom pokrenuti sve cells. *)

Modules

DBSCAN algorithm

Input: data set 𝒜 ⊂ (ℝ^n), MinPts, ϵ-density ;The algorithm can use data sets from E-ColiA.txt;
Output: partition Π of the set 𝒜 with marked border points SSE and outliers SSN;
Submodules: Neps[], EPSILON[], CorePoint[], DBSCAN1[];

DBSCAN[A_,MinPts_,eps_,Ind_:0]:=Module[{kk=0,It0=.1Length[A],TH=.1 Length[A],time=0,B,cp,PPP,PP,SPP,SE,SSN},
PPP=SSN=SSE=SPP={};
B=A; TH=.01 Length[A];
While[Length[B]>TH MinPts,kk=kk+1;cp=CorePoint[B,eps,It0];If[Ind≠0,Print["cp: ",cp]];  
dbran=DBSCAN1[B,cp[[1]],MinPts,eps]//Timing;
PP=dbran[[2,1]];Border=dbran[[2,2]]; 
SPP=Append[SPP,PP];
Noise=dbran[[2,3]];  
time=time+dbran[[1]]; 
B=Complement[B,Union[PP,Border,Noise]]; 
If[Ind≠0,Print[kk,". izabrana točka cp: ",cp[[1]],"; |Cp|=",cp[[2]],"; |PP|=",Length[PP],"; |Border|= ",Length[Border],"; |B|=",Length[B],"; Time(DBRAN)=",time]];
 (* Izabrani cluster + Ostatak skupa - može se uključiti 
slizbor=ListPlot[B,PlotStyle->{Blue,PointSize[.005]},AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,PlotLabel->"Ostatak"];
slPP=ListPlot[PP,PlotStyle->{Blue,PointSize[.001]},AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,PlotLabel->"Izabrani klaster"];
Print[GraphicsGrid[{{slPP,slizbor}},ImageSize->Medium]];
*)
If[Ind≠0,Print[Show[ListPlot[Flatten[SPP,1],PlotStyle->{Gray,PointSize[.001]}],ListPlot[B,PlotStyle->{Blue,PointSize[.005]}],AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG ,Frame->True,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]]];
PPP=AppendTo[PPP,PP]; SSE=Append[SSE,Border];SSN=Append[SSN,Noise];
];
SSN=Complement[A,Flatten[PPP,1]];
Print["|PPP| = ",Sum[Length[PPP[[j]]],{j,Length[PPP]}]," od toga |SSE| = ",Sum[Length[SSE[[j]]],{j,Length[SSE]}],";  |SSN| = ",Length[SSN], ";  TOTAL[A] = ",Length[A],"; Total[time]=",time];
slPPP=ListPlot[Flatten[PPP,1]];
slSN=ListPlot[SSN,PlotStyle->{Red,PointSize[.025]}];
slSE=ListPlot[Flatten[SSE,1],PlotStyle->{Brown,PointSize[.02]}];
sl1=ListPlot[PPP, PlotStyle->{PointSize[.015]},AspectRatio->Automatic,Frame->True,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[sl1];
sl2=Show[sl1,slSN,slSE,AspectRatio->Automatic,Frame->True,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
Print[GraphicsGrid[{{sl1,sl2}},ImageSize->Medium]];
{PPP,SSE,SSN}
]

A single step of DBSCAN algorithm

Input: data set 𝒜 ⊂ (ℝ^n), core point (Subscript[c, p]), MinPts, ϵ-density ϵ(𝒜);
Output: cluster π, border points SSE, outliers SSN;
Submodules: Neps[], EPSILON[]

DBSCAN1[A_,xx_,MinPts_,eps_]:=Module[{B,S,SE,SN,PP},
SE={};
B=A;
S={xx}; SN={};PP={}; 
    Ne=Neps[A,xx,eps];
If[Length[Ne]>1&&Length[Ne]< MinPts,SE=Union[SE,Ne]];
    If[Length[Ne]< MinPts,SN=Union[SN,{xx}],
         S=Union[S,{xx}];
	NeB=Intersection[B,Ne];
        If[NeB≠{},PP=Union[PP,NeB];
	B=Complement[B,NeB]];
	Qu=Complement[Ne,S];
         While[Qu≠{},y=RandomChoice[Qu]; S=Union[S,{y}];
                   Ny=Neps[A,y,eps];
If[Length[Ny]>1&&Length[Ny]< MinPts,SE=Union[SE,Ny]];
                   If[Length[Ny]≥ MinPts,NeB=Intersection[B,Ny];
                        PP=Union[PP,NeB];
		      B=Complement[B,Ny];      
                        SN=Complement[SN,Ny]; 
                        Qu=Union[Qu,Complement[Ny,S]];
                      ];
                    Qu=Complement[Qu,{y}];
                 ];  
      ];
{PP,SE,SN,xx}
]

ϵ-density ϵ(𝒜) of the set 𝒜 ⊂ (ℝ^n)

Input: data set 𝒜 ⊂ (ℝ^n), MinPts ≥ 2;
Output: ϵ-density ϵ(𝒜);
Submodules: dM[]

EPSILON[A_,MinPts_,Fig_:0]:=Module[{yt,yy,ee,eps,slline,sleps},
yt=Table[
{A[[i]],Max[Table[Norm[Nearest[A,A[[i]],MinPts][[j]]-A[[i]]],{j,MinPts}]//N]}
,{i,Length[A]}]; 
yy=yt[[All,2]]; pp=95.;
eps=Quantile[yy,pp/100];
If[Fig≠0,ee={Quantile[yy,50/100],Quantile[yy,92/100],eps,Quantile[yy,99.5/100]};
slline=Plot[ee,{xx,0,Length[yy]},PlotStyle->{Black,Brown, Red,Orange},PlotRange->All];
sleps=Print[Show[ListPlot[yy,PlotStyle->{PointSize[.008]},PlotRange->All],slline,ImageSize->Small]]];
eps
]

Set (Subscript[N, 𝒜])(p,ϵ)={q ϵ 𝒜 : d(q,p) < ϵ} = B(p,ϵ)⋂𝒜 where B(p,ϵ)⋂𝒜 ⊂  (ℝ^n)is the open ϵ-ball around p

Input: data set 𝒜 ⊂ (ℝ^n), p ϵ 𝒜, ϵ > 0;
Output: (Subscript[N, 𝒜])(p,ϵ);
Submodules: none needed

Neps[A_,x_,eps_]:=Nearest[A,x,{All,eps},DistanceFunction->(Norm[#1-#2,2]&)]

Choice of the core point in the set 𝒜 , based on  (Subscript[It, 0]) random choices

Input: data set 𝒜 ⊂ (ℝ^n), ϵ > 0, number of random choices (Subscript[It, 0]);
Output: core point (Subscript[c, p]), number of points |(Subscript[C, p])| in ϵ-neighborhood of the core point (Subscript[c, p]);
Submodules: none needed

CorePoint[A_,eps_,It0_]:=Module[{kk=1,rc,Ne,B,br,br1,rep,Rep},
rc=RandomChoice[A]; Ne=Neps[A,rc,eps]; br=Length[Ne];
(* rep - chosen point; Rep -  number of points in the neighborhood of point rep *)
rep=rc;Rep=Ne;kk=1;B=Complement[A,{rc}];
While[kk<It0,kk=kk+1;
rc=RandomChoice[B];Ne=Neps[B,rc,eps]; br1=Length[Ne];
If[br1>br, br=br1;rep=rc;Rep=Ne]
	];
{rep,Length[Rep]}
]

Distance from a point to the segment in the plane

(* Input: point T ∈ (ℝ^2), endpoints μ,ν of the segment [μ,ν]; 
Output: distance;
Submodules: none needed *)

dist[t_,Seg_]:=Module[{lam,a=Seg[[1]],b=Seg[[2]],t0,dd},
If[(b-a).(t-a)<0,t0=a;dd=Norm[t-a],
If[(a-b).(t-b)<0,t0=b;dd=Norm[t-b],
lam=(b-a).(b-t)/Norm[b-a]^2;
t0=lam a+(1-lam)b;
dd=Norm[t0-t]]];  
dd
];

Dmn[t_,xm_,ym_,xn_,yn_]:=Module[{lam,a,b,t0,dd},
a={xm,ym};b={xn,yn};
If[(b-a).(t-a)<0,t0=a;dd=Norm[t-a],
If[(a-b).(t-b)<0,t0=b;dd=Norm[t-b],
lam=(b-a).(b-t)/Norm[b-a]^2;
t0=lam a+(1-lam)b;
Norm[t0-t]]]
];

Udaljenost točke T do G - circle sa središtem u Seg radijusa r

(* Seg je segmen (mu,nu) *)
DGcir1[T_,Seg_,r_]:=Module[{},Return[Abs[dist[T,Seg]-r]]
                                             ];
DGcir[T_,Seg_,r_]:=Module[{},
(dist[T,Seg]^2-r^2)^2
                                             ];

Objective function for Generalized circles

WGLS[PI_,Circ_]:=Module[{},
Return[
Sum[Sum[PI[[j,s,1]]DGcir[PI[[j,s,2]],Circ[[j,1]],Circ[[j,2]]],{s,Length[PI[[j]]]}]
,{j,Length[PI]}]
]]

Weighted k - closest generalized circles algorithm (WKGC)

(* Circ su IA za GC 
Koristi se DIRECT algorithm without Compile *)
WKGC[w_,A_,Circ_,Iter_]:=Module[{It=0,m=Length[w],k=Length[Circ],Kr,PI,tab,fc,min,imin,Pod,F0,F1,p,q,r,dom,min1,fCompile,FF,EPS=0.5},
Kr=Table[0,{j,k}];
Pod=Table[{w[[i]],A[[i]]},{i,m}];
						(* Odredjivanje particije *)
PI=Table[{},{i,k}];
Do[
tab=Table[DGcir[Pod[[i,2]],Circ[[j,1]],Circ[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F1=WGLS[PI,Circ]; 
                                              (* Računanje kružnica-centara *)
Do[
c0=Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]];
Cov=Sum[PI[[j,s,1]]KroneckerProduct[(PI[[j,s,2]]-c0),(PI[[j,s,2]]-c0)],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
v=Eigensystem[Cov][[2,1]];
d0=2Sum[Abs[Dot[PI[[j,s,2]]-c0,v]],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
seg0={c0-.5d0 v,c0+.5d0  v};
r0=Mean[Table[dist[PI[[j,s,2]],seg0],{s,Length[PI[[j]]]}]];
(* Print[j,". appr: ",{seg0,r0}]; *)
lam=d0/4;
dom={{seg0[[1,1]]-lam,seg0[[1,1]]+lam},{seg0[[1,2]]-lam,seg0[[1,2]]+lam},{seg0[[2,1]]-lam,seg0[[2,1]]+lam},{seg0[[2,2]]-lam,seg0[[2,2]]+lam},{r0-lam,r0+lam}};(* Print[j,". dom: ",dom]; *)
distGC[xm_,ym_,xn_,yn_,r_]:=Sum[(Dmn[PI[[j,s,2]],xm,ym,xn,yn]^2-r^2)^2,{s,Length[PI[[j]]]}];
(*FCompile=Compile[{{z1,_Real},{z2,_Real},{z3,_Real},{z4,_Real},{z5,_Real}},
Evaluate[Sum[(Dmn[PI[[j,s,2]],z1,z2,z3,z4]^2-z5^2)^2,{s,Length[PI[[j]]]}]]];	
distGC[z1_,z2_,z3_,z4_,z5_]:=Apply[FCompile,{z1,z2,z3,z4,z5}];*)
(*Print[distGC[j,seg0[[1,1]],seg0[[1,2]],seg0[[2,1]],seg0[[2,2]],r0]];*)
tDIR=Timing[min=DIRECT[distGC,dom,DMin->10.^(-4),MaxIter->30, Eps->10^(-4)]];
Kr[[j]]={{{min[[1,1]],min[[1,2]]},{min[[1,3]],min[[1,4]]}},min[[1,5]]};
,{j,k}];
F0=WGLS[PI,Kr];  
While[
(* Chop[F0-F1]≠0&&It<Iter, *)
Abs[F0-F1]/F0>EPS&&It<Iter,
F0=F1;It=It+1;
PI=Table[{},{i,k}];
Do[
tab=Table[DGcir[Pod[[i,2]],Kr[[j,1]],Kr[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
Kr=Table[0,{j,k}];
Do[
c0=Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]];
Cov=Sum[PI[[j,s,1]]KroneckerProduct[(PI[[j,s,2]]-c0),(PI[[j,s,2]]-c0)],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
v=Eigensystem[Cov][[2,1]];
d0=2Sum[Abs[Dot[PI[[j,s,2]]-c0,v]],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
seg0={c0-.5d0 v,c0+.5d0  v};
r0=Mean[Table[dist[PI[[j,s,2]],seg0],{s,Length[PI[[j]]]}]];
(* Print[j,". appr: ",{seg0,r0}]; *)
lam=d0/4;
dom={{seg0[[1,1]]-lam,seg0[[1,1]]+lam},{seg0[[1,2]]-lam,seg0[[1,2]]+lam},{seg0[[2,1]]-lam,seg0[[2,1]]+lam},{seg0[[2,2]]-lam,seg0[[2,2]]+lam},{r0-lam,r0+lam}};(* Print[j,". dom: ",dom]; *)
distGC[xm_,ym_,xn_,yn_,r_]:=Sum[(Dmn[PI[[j,s,2]],xm,ym,xn,yn]^2-r^2)^2,{s,Length[PI[[j]]]}];
(*FCompile=Compile[{{z1,_Real},{z2,_Real},{z3,_Real},{z4,_Real},{z5,_Real}},
Evaluate[Sum[(Dmn[PI[[j,s,2]],z1,z2,z3,z4]^2-z5^2)^2,{s,Length[PI[[j]]]}]]];	
distGCz1_,z2_,z3_,z4_,z5_]:=Apply[FCompile,{z1,z2,z3,z4,z5}];*)
(*Print[distGC[j,seg0[[1,1]],seg0[[1,2]],seg0[[2,1]],seg0[[2,2]],r0]];*)
tDIR=Timing[min=DIRECT[distGC,dom,DMin->10.^(-4),MaxIter->15, Eps->10^(-4)]];
Kr[[j]]={{{min[[1,1]],min[[1,2]]},{min[[1,3]],min[[1,4]]}},min[[1,5]]};
(* Print[Kr[[j]]];*)
,{j,k}];
F1=WGLS[PI,Kr];  
];
{PI,Kr,F1,It}
]

Density Based Clustering - Index

DBC[PI_,Circles_]:=Module[{k1,h,rho={},mad={},QAD={},rj},
k1=Length[Circles]; 
Do[
mj=Length[PI[[j]]];
dC=2Pi Circles[[j,2]]+2Norm[Circles[[j,1,1]]-Circles[[j,1,2]]];
ud1=Table[DGcir1[PI[[j,s,2]],Circles[[j,1]],Circles[[j,2]]],{s,mj}];
QAD=Append[QAD,Mean[ud1]+2Sqrt[Variance[ud1]]];
,{j,Length[Circles]}]; 
Print["Veličine ",k1," klastera: ",LPI=Table[Length[PI[[j]]],{j,k1}]];
(*Print["QAD: ",QAD, "; Max[QAD] = ",Max[QAD]]; *)
{Max[QAD],QAD}
];

DIRECT algorithm for global optimization

Input: objective function f:dom⟶ℝ, dom ⊂ (ℝ^n), domain dom of
the function f, additional options: DMin⟶(10^(-3))(minimal length of rectangle sides), MaxIter⟶50 (maximal number of iteration), Eps⟶(10^(-4))(parameter ϵ from the method);
Output:  {(x^⋆),f((x^⋆))}, where (x^⋆)∈ (Underscript[argmin, x∈ dom]) f(x);
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,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];
]

Data

figureA={,,,};
SetDirectory[NotebookDirectory[]]; 
SlikeA=ReadList["E-ColiA.txt"];
Length[SlikeA]

ii=1; figure=figureA[[ii]];
A=SlikeA[[ii]]//N;  
m=Length[A];
w=Table[1,{i,m}];
del1=del2=10;
xG={Min[A[[All,1]]]-del1,Max[A[[All,1]]]+del1};
yG={Min[A[[All,2]]]-del2,Max[A[[All,2]]]+del2};
AG={xG,yG};
MinPts=Floor[Log[m]];  

eps=EPSILON[A,MinPts,0];
Print["|A|=",m,";  MinPts= ",MinPts,";  eps=",eps,";  Lower Bound(rho)=",LB=MinPts/(2 eps),"\nDimenzija slike: ",AG];
slA=slA0=ListPlot[A,PlotStyle->{PointSize[.02]},AspectRatio->Automatic,Frame->True,GridLines->Automatic,PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]
p=2;  TINI=TDIR={};

Searching for an initial partition - call DBSCAN

db=DBSCAN[A,MinPts,eps,0];

Call WKGC

(* Traženje initial approximation za GC *)
PIP=db[[1]]; k=Length[PIP];
PI=Kr=Kr0=Table[0,{j,k}];TINI=TDIR={};
Do[
PI[[j]]=Table[{1,PIP[[j,s]]},{s,Length[PIP[[j]]]}];
,{j,k}]
Do[
tINI=Timing[
c0=Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]];
Cov=Sum[PI[[j,s,1]]KroneckerProduct[(PI[[j,s,2]]-c0),(PI[[j,s,2]]-c0)],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
v=Eigensystem[Cov];
d0=2Sum[Abs[Dot[PI[[j,s,2]]-c0,v[[2,1]]]],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
seg0={c0-.5d0 v[[2,1]],c0+.5d0 v[[2,1]]};
(* r0 je obicni Mean *)
r0=Mean[Table[dist[PI[[j,s,2]],seg0],{s,Length[PI[[j]]]}]];
(* r0 prema Thomas2011 *)
r0=Sum[Abs[Dot[PI[[j,s,2]]-c0,v[[2,2]]]],{s,Length[PI[[j]]]}]/Length[PI[[j]]];
];
Print[Style["Initial approximation for j=",14,Blue],Style[j,14,Blue],": ",Kr0[[j]]={seg0,r0}];
(* Application of DIRECT algorithm *)
lam=d0/4;
dom={{seg0[[1,1]]-lam,seg0[[1,1]]+lam},{seg0[[1,2]]-lam,seg0[[1,2]]+lam},{seg0[[2,1]]-lam,seg0[[2,1]]+lam},{seg0[[2,2]]-lam,seg0[[2,2]]+lam},{r0-lam,r0+lam}};(* Print[j,". dom: ",dom]; *)
Clear[xm,ym,xn,yn,r];
distGC[xm_,ym_,xn_,yn_,r_]:=Sum[(Dmn[PI[[j,s,2]],xm,ym,xn,yn]^2-r^2)^2,{s,Length[PI[[j]]]}]; 
(*mj=Length[PI[[j]]];
FCompile=Compile[{{z1,_Real},{z2,_Real},{z3,_Real},{z4,_Real},{z5,_Real}},
Evaluate[Sum[(Dmn[PI[[j,s,2]],z1,z2,z3,z4]^2-z5^2)^2,{s,mj}]]];	
distGC[z1_,z2_,z3_,z4_,z5_]:=Apply[FCompile,{z1,z2,z3,z4,z5}]; *)
(*Print[distGC[j,seg0[[1,1]],seg0[[1,2]],seg0[[2,1]],seg0[[2,2]],r0]];*)
tDIR=Timing[min=DIRECT[distGC,dom,DMin->10.^(-4),MaxIter->25, Eps->10^(-4)]];
Kr[[j]]={{{min[[1,1]],min[[1,2]]},{min[[1,3]],min[[1,4]]}},min[[1,5]]};
Print[Style["Solution by DIRECT for j=",14,Blue],Style[j,14,Blue],": ",Kr[[j]],"\n"];
TINI=Append[TINI,tINI[[1]]];TDIR=Append[TDIR,tDIR[[1]]];
,{j,k}]
					(* FIGURE IA *)
Print[Style["Figure Initial approximation:",14,Blue]];
slSeg=Table[Graphics[{Red,Thick,Line[{Kr0[[j,1]]}]}],{j,Length[Kr0]}];
slGcir=Table[ContourPlot[dist[{x,y},Kr0[[j,1]]]==Kr0[[j,2]],{x,AG[[1,1]],AG[[1,2]]},{y,AG[[2,1]],AG[[2,2]]},ContourStyle->{Red},AspectRatio->Automatic],{j,Length[Kr0]}];
Print[Show[sl1,slGcir,slSeg,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]];
					(* FIGURE DIRECT *)
Print[Style["Figure DIRECT:",14,Blue]];
slSeg=Table[Graphics[{Red,Thick,Line[{Kr[[j,1]]}]}],{j,Length[Kr]}];
slGcir=Table[ContourPlot[dist[{x,y},Kr[[j,1]]]==Kr[[j,2]],{x,AG[[1,1]],AG[[1,2]]},{y,AG[[2,1]],AG[[2,2]]},ContourStyle->{Red},AspectRatio->Automatic],{j,Length[Kr]}];
Print[Show[sl1,slSeg,slGcir,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]];
								(* WKGC *)
tWGKC=Timing[sol=WKGC[w,A,Kr,50]];
						(* Results and Figures *)
Print[Style["Figure WKGC:",14,Blue]];
Print[Style["Solution by WKGC: ",14,Blue],sol[[2]]];
slSeg=Table[Graphics[{Red,Thick,Line[{sol[[2,j,1]]}]}],{j,Length[sol[[2]]]}];
slGcir=Table[ContourPlot[dist[{x,y},sol[[2,j,1]]]==sol[[2,j,2]],{x,AG[[1,1]],AG[[1,2]]},{y,AG[[2,1]],AG[[2,2]]},ContourStyle->{Red},AspectRatio->Automatic],{j,Length[sol[[2]]]}];
Print[detection=Show[sl1,slGcir,slSeg,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]];
Print[Style["QAD and QUAD per generalized circles: ",14,Blue],Style[DBC[sol[[1]],sol[[2]]],14,Blue]];
Print[Style["Total time Non-Compiled {INI,DIR,WGKC}: ",14,Blue],Style[{TINI,TDIR,tWGKC[[1]]},14,Blue]];
Print[Style["CPU times for Compiled evaluations is obtained dividing with approximately 150",12,Orange]];

(* za A: del1=0; del2=48;
za B: del1=1; del2=38;
*)

Comparison

del1=0;del2=48;cen=Table[sol[[2,j,1]],{j,Length[sol[[2]]]}]-Table[{{del1,del2},{del1,del2}},{j,Length[sol[[2]]]}]
slSeg=Table[Graphics[{Red,Thick,Line[cen[[j]]]}],{j,Length[sol[[2]]]}];
slGcir=Table[ContourPlot[dist[{x,y},cen[[j]]]==sol[[2,j,2]],{x,0,130},{y,0,150},ContourStyle->{Red},AspectRatio->Automatic],{j,Length[sol[[2]]]}];
detection=Show[slGcir,slSeg,AspectRatio->Automatic,GridLines->Automatic,AxesOrigin->{0,0},ImageSize->Small];
Show[figure,detection,ImageSize->Small]