Trennung von Punktmengen mit Linearer Optimierung (nicht linear)

Funktion nimmt Punktmengen A und B entgegen; zusätzlich maximale Schritte als Steps; OutPut sind die einzelnen Trennungsebenen {{c,alpha,beta}}. Der erste Eintrag ist dabei die Ebenen nach dem 1. Schritt.

LP[AAA_,BBB_,steps_] := Module[{m,k,A,B,F,DD,d,b,x,cc,L,n,AA,BB,aT,bT,cT,j,alpha, beta,II,JJ,i,l,ifert,EE,I1,imax},
imax=steps;
A=AAA;
B=BBB;
n=Length[A[[1]]]; (* Dimension *)
F=IdentityMatrix[n]; (* Einheitsmatrix Dimension n *)
p=0;
For[i=1,i<=imax,i++,p=p+1;
I1[x_]:= If[x==0,Table[1,2 n],ReplacePart[Table[1,2 n],x->-1]] ;(* Vektor mit 2 n mal 1 für b *)
m=Length[A]; (* Anzahl Muster A *)
k=Length[B]; (* Anzahl Muster B *)
DD=ArrayFlatten[{{-A,1,0,A,-1,0},{B,0,-1,-B,0,1},{F,0,0,-F,0,0},{-F,0,0,F,0,0}}]; (* Bedingungsmatirx *)
d=Join[Table[0,n],{1},{-1},Table[0,n],{-1},{1}]; (* Zielfunktionsvektor *)
b=Join[Table[0,m],Table[0,k],I1[0]];
x=LinearProgramming[-d,-DD,-b];
cc=Take[x,n+2]-Take[x,-(n+2)];
If[x==Table[0,2 n +4], (* keine lineare Trennung *)
L=10^9;
AA={};
BB={};
aT[i]=0;
bT[i]=0;
cT[i] = Table[0,n]; 
(* keine lineare Trennung *)
For[j=1,j<=2 n, j++, 
b=Join[Table[0,m],Table[0,k],I1[j]]; (* Restriktionsvektor *)
x=LinearProgramming[-d,-DD,-b];(* Lösung von min -d.x, -DD.x >= -b, x>= 0 = max d.x, DD.x <= b, x >=0 *)
cc=Take[x,n+2]-Take[x,-(n+2)]; (* Ermittlung von c,alpha, beta *)
c = Take[cc,n];
alpha = Take[cc,-2][[1]];
beta = Take[cc,-1][[1]];
II={};JJ={};
For[l = 1, l <= m, l++, If[c . A[[l]] - beta < 0, II = Append[II, A[[l]]]]];
For[l = 1, l <= k, l++, If[c . B[[l]] - alpha > 0, JJ = Append[JJ, B[[l]]]]];
If[Length[II]+Length[JJ]<L,
L=Length[II]+Length[JJ];AA=II;BB=JJ;
aT[i]=alpha;bT[i]=beta;cT[i]=c;
]];A=AA;B=BB,
cT[i] = Take[cc,n];
aT[i] = Take[cc,-2][[1]];
bT[i] = Take[cc,-1][[1]];
i=imax;];
EE={};
EE=Table[{cT[l],aT[l],bT[l]},{l,1,p}]];Return[EE]];

TestLP[P,EE,n] ist eine Testfunktion für unbekannte Punkte; Input ist der Punkt P und die Ebenen EE={{c,alpha,beta}} (Output aus der Funktion LP[A,B,step]). n gibt an, bis welcher Ebene geteset werden soll (max. n=Length[EE]).  Output ist die Zuordnung des Punktes P auf A bzw. B und auf welcher Ebene der Punkt gefunden wurde (Output z.B. {A,2}). Je höher die Ebene, umso unsicherer kann der Punkt P zugeordnet werden.

TestLP[P_,EE_,n_]:=Module[{i,K,k,t},
t=10^-9;
For[i=1,i<=n,i++,
If[EE[[i]][[1]].P-EE[[i]][[3]]-t>=0,K="A",If[EE[[i]][[1]].P-EE[[i]][[2]]+t<=0,K="B", K="N";k=0]];
If[K!="N",k=i;i=n]];
Return[{K,k}]]

RestLP[AA,BB,EE,k] gibt an, welche Punkte nach k Trennungen (Schritte) des Algortihmus LP[AA,BB,step] noch übrig sind. Output ist {RestA, RestB}.

RestLP[AA_,BB_,EE_,k_]:= Module[{a,b,AAA,BBB,i,n},
AAA={};
BBB ={};
n=Length[EE];
For[i=1,i<=Length[AA],i++,
If[TestLP[AA[[i]],EE,n][[2]]>=k+1,AAA=Append[AAA,AA[[i]]]]];
For[i=1,i<=Length[BB],i++,
If[TestLP[BB[[i]],EE,n][[2]]>=k+1,BBB=Append[BBB,BB[[i]]]]];
Return[{AAA,BBB}]]

Plot der Trennung (2D) für Punkte A und Punkte B mit der Trennung cc={{c},alpha,beta}. Falls f=1 wird auch die Ebene zu (alpha+beta)/2 geplotet. g gibt die Größe der gezeichneten Punkte an (z.B. g=0.05)

PlotLP2D[A_,B_,cc_,f_,g_]:= Module[{m,k,alpha,beta,mAx,mIx,mIy,c,cH,d,e,aa,bb,AA1,dd},
m=Length[A];
k=Length[B];
alpha=cc[[2]];
beta =cc[[3]];
c=cc[[1]];
mAx=Max[Max[Transpose[A][[1]]], Max[Transpose[B][[1]]]] (* Maximum x in Punktmenge *);
mIx=Min[Min[Transpose[A][[1]]], Min[Transpose[B][[1]]]] (* Minimum x in Punktmenge *);
mAy=Max[Max[Transpose[A][[2]]], Max[Transpose[B][[2]]]](* Maximum y in Punktmenge *);
mIy=Min[Min[Transpose[A][[2]]], Min[Transpose[B][[2]]]](* Minimum y in Punktmenge *);
cH=Plot[((alpha+beta)/2 -c[[1]]*x)/c[[2]],{x,mIx,mAx},AxesStyle-> {Directive[Red,26]},PlotStyle->{RGBColor[140/255,35/255,0]}]; (* Hyperebene *)
d=Plot[(alpha -c[[1]]*x)/c[[2]],{x,mIx,mAx},PlotStyle->{RGBColor[140/255,35/255,0],Dotted},AxesStyle-> {Directive[Thick,26]}]; (* Supportebene Punktmenge A *)
e=Plot[(beta -c[[1]]*x)/c[[2]],{x,mIx,mAx},PlotStyle->{RGBColor[140/255,35/255,0],Dashed},AxesStyle-> {Directive[Thick,26]}]; (* Supportebene Punktmenge B *)
aa=Table[ Graphics[Circle[A[[a]],g]],{a,m}]; (* Punkmenge A *)
bb=Table[Graphics[Disk[B[[a]],g]],{a,k}];(* Punkmenge B *)
AA1=Show[{aa,bb,If[f==1,cH,{}],d,e},PlotRange->{mIy,mAy},Axes->True,LabelStyle-> {FontFamily -> "Courier"}, AxesStyle-> {Directive[Thickness[0.006],14],Directive[Thickness[0.006],14]}];
Return[AA1]]

Beispiel:

A={{-5.853944095`,20.50080539`,-2.922694725`,6.577985675`,-7.304068366`,16.08069586`,-9.486973655`,-7.304601899`,-3.65497076`},{-20.06880734`,-10.77276645`,7.908548422`,1.436781609`,-5.746300819`,-8.624407072`,3.596604805`,-0.719062343`,5.03379836`},{-5.790807094`,10.86405447`,-13.02271741`,-24.63054187`,-2.90486565`,-4.358564579`,13.80813953`,2.177226214`,-5.804672762`},{1.462308986`,8.041523503`,5.843681519`,12.41057089`,2.916514765`,5.831328814`,-9.470386829`,-1.458363716`,-1.458576429`},{-13.19841619`,-11.01321586`,-11.76038221`,2.207667967`,4.414361389`,-7.35402265`,13.24698263`,0.734969866`,0},{0.713165026`,16.4016259`,-1.423892923`,-2.136143549`,15.66839969`,0,-5.688686624`,-6.403415155`,7.119464616`},{10.48610591`,4.489337823`,11.21830828`,2.988196623`,-1.493651979`,4.481625336`,-7.466029565`,-2.988643156`,-0.747384155`},{-8.444757213`,-9.156219186`,-4.229820233`,1.410536709`,2.115506664`,-7.755217146`,20.46144077`,3.520630897`,-2.815513479`},{-9.107468124`,-3.798237618`,5.319553158`,-6.076257026`,-2.27998176`,-8.361839605`,9.890444309`,2.280155051`,7.598784195`},{-8.278145695`,-0.753182195`,-3.766194637`,0.753522719`,2.26039783`,-6.026365348`,6.029999246`,7.532956685`,4.516371848`}};

B={{-2.224529141`,-2.225024104`,-1.483679525`,5.193648909`,0.741564702`,-2.966038855`,-6.675567423`,3.711125956`,-2.967799377`},{3.675389591`,-1.469615695`,8.084074373`,-6.608899985`,1.469615695`,3.673499376`,-3.672150411`,3.673499376`,1.468860165`},{0.735186002`,2.940527825`,2.204747556`,2.204261572`,7.345919342`,-5.138368935`,0.734430082`,0.734376147`,2.202966662`},{0,1.471345546`,-2.942258183`,2.207343095`,6.620567898`,-2.205395869`,-3.676470588`,-1.471129092`,6.621054955`},{-1.46113384`,5.114715768`,1.460600307`,2.190580504`,-7.300335815`,-6.575102279`,4.386285547`,-10.9609061`,-1.463057791`},{5.8317539`,2.185633105`,-2.18515551`,-5.099810578`,2.186748305`,-3.643783705`,1.458044762`,7.289161018`,-6.555466531`},{-16.62089897`,-6.51465798`,-7.243227582`,2.899391128`,2.173913043`,9.418242411`,-2.895193978`,1.448016218`,2.895613146`},{-13.74222479`,2.172810893`,4.34467777`,0,-0.723798494`,-7.962359754`,-9.417560127`,-2.900442317`,7.253209545`},{-11.51907847`,3.603863342`,-4.323078032`,6.48742161`,-9.364644864`,-3.605162593`,2.163877669`,-5.769092089`,-3.607763908`},{1.42969476`,-4.28847116`,8.580622095`,-4.28663285`,3.573725967`,5.00142898`,2.856530743`,-9.281073749`,-18.57939117`}};

EE=LP[A,B,7]

TestLP[B[[7]],EE,2]

EE[[1]][[1]].A[[6]]

EE[[1]][[1]]

{-0.08789971115265688`,0.32690829601483906`,0.4687061816636993`,0.22899838798230343`,0.07329250385805898`,0.6566117754079289`,1.`,0.08556857328876703`,0.32355668024323825`}.{0.713165026`,16.4016259`,-1.423892923`,-2.136143549`,15.66839969`,0,-5.688686624`,-6.403415155`,7.119464616`}

A[[6]]

Table[TestLP[B[[i]],EE,1],{i,1,Length[B]}] (* nach ersten Schritt klassifizierte A-Punkte *)

Table[TestLP[A[[i]],EE,1],{i,1,Length[A]}] (* nach zweitem Schritt klassifizierte A-Punkte *)

A1=RestLP[A,B,EE,1][[1]] (* Rest der Punkte A nach Schritt 1 *)
B1 =RestLP[A,B,EE,1][[2]] (* Rest der Punkte B nach Schritt 1 *)

AA=PlotLP2D[A,B,EE[[1]],0,0.05]

Export["C:\\MyTex\\LaTeXTemplates_wenneker-article_v2.0\\Trennung2.eps",AA]