Numerické hledání Steinerova bodu pro konvoluci s kruhovým jádrem

V tomto souboru se nachází definice potřebných funkcí (definována je množina A, konvoluční jádro, funkce příslušnosti fuzzy množiny vzniklé konvolucí s daným jádrem apod.). 

Kromě definic jsou tu funkce pro výpočet bodů na alpha hladině i numerickou integraci pro výpočet Steinerova bodu. 

Na konci je rozsáhlá vizualizace. 

(Před voláním funkce doItAll je třeba inicializovat všechny funkce.)
(soubory jsou ukladany / nacitany do / z domovskeho adresare - resp. do / z defaultne nastavene cesty v programu Mathematica)

Nastavení

Množina A

(* geometricky objekt, pres ktery se bude integrovat -- volano jako Element[{x,y},setA] *)

setA = Disk[{0,0},1,{0,Pi}];

Konvoluce s vybraným jádrem

Nastavení parametrů jádra

R = 1/2;
PARDSC = "R1_2";
(* kernel[x,y] je funkce jadra pouzita pri konvoluci *)
kernel[x_?NumericQ,y_?NumericQ] = Piecewise[{{1/(Pi (R^2)),(x^2)+(y^2)≤ (R^2)}}];

Konvoluce a její maximum

(* konvoluce z definice jako prunik*)
convolution[x_ ?NumericQ,y_?NumericQ]:=intersect[x,y,R,0,0]; (* alternativne regIntersect[x,y,diskR,A] *)
(* maximum -- prevracena hodnota pouzita jako normalizacni clen *)
{muStartX,muStartY} = {0,1/2};
maximum = convolution[muStartX,muStartY];

Funkce příslušnosti fuzzy množiny vzniklé konvolucí

fuzzyMemFunction[x_?NumericQ,y_?NumericQ] := (1/maximum) * convolution[x,y];
fuzzyMemFunction::usage = "fuzzyMemFunction[x_?NumericQ,y_?NumericQ] vrati stupen prislusnosti daneho bodu (x,y) k nasi fuzzy mnozine";

Další parametry použité v průběhu výpočtu (globální proměnné):

kernelName = "Disk"; (* nazev jadra - pro pojmenovano souboru *)
minPoints = 32; (* minimalni pocet bodu hledanych na α-hladine *)
maxPoints = 128; (* maximalni pocet bodu hledanych na α-hladine *)
PHI = 0; (* pootoceni pocatecniho smeru pri generovani smerovych vektoru, volano v generateCirclePoints *)
PRECISION = 100; (* presnost pri volani N[..., PRECISION] *)
NUMBERofPOINTS = maxPoints; (* pocet bodu hledanych na α-hladine *)
EPSILON = (10^(-10)); (* maximalni chyba pri hledani  bodu na hladinach *)
maxITER = 500; (* maximalni pocet iteraci pri hledani 1 bodu na α-hladine *)
integrationMETHOD = NrombergY;(* vyber integracni metody *)
(* NrombergX[NoScale] NrombergY[NoScale] NsimpsonY NtrapezoidX[NoScale] NtrapezoidY[NoScale] *)
integrationMETHODNS = NrombergYNoScale;(* poze NoScale metody*)
alphaSCALE = 1.5; (* skalovaci faktor pro alpha (krok) pri hledani bodu na hladine *)

Nastavení jmen souborů a hladin, které se budou hledat

(* hledane vysky α *)
 HEIGHTS = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9};

(* popisky souboru pro body na hladinach *)

heightsLABELS = Map[StringJoin["LAYER-",kernelName,"-",PARDSC,"-",ToString[NUMBERofPOINTS],
"-",ToString[#,InputForm],".dat"]&,HEIGHTS] // StringReplace[#,"/":>"_"]&;
(* popisky souboru pro Steinera *)
steinerLABELS = 
Map[StringJoin["STEINER-",kernelName,"-",PARDSC,"-",ToString[NUMBERofPOINTS],
"-",ToString[#,InputForm],".dat"]&,HEIGHTS] // StringReplace[#,"/":>"_"]&;

Základní vizualizace fuzzy množiny

Pokud byly funkce v sekci nastavení inicializovány správně, následuje obrázek fuzzy množiny a jejích vybraných hladin:

Fuzzy množina

Plot3D[regIntersect[x,y,R,setA],{x,-2,2},{y,-2,2},PlotRange->All]

Vybrané hladdiny

ContourPlot[regIntersect[x,y,R,setA],{x,-2,2},{y,-2,2},Contours->{0.2,0.4,0.6,0.8},AspectRatio->Automatic,PlotRange->All]

Výpočet

doItAll[HEIGHTS,heightsLABELS,steinerLABELS,NUMBERofPOINTS, EPSILON,maxITER]

Funkce nutné pro výpočet

Generování směrových vektorů

generateCirclePoints[number_] := Table[{Cos[x + PHI],Sin[x + PHI]},{x,0,2Pi,2*Pi/number}][[;;-2]];
generateCirclePoints::usage="generateCirclePoints[number_] trati tabulku [number] rovnomerne rozdelenych souradnic smerovych vektoru";

Funkce pro výpočet plochy průniku

Průsečíky kružnice Circ[{u,v},radius] s úsečkou [{shiftX-1,shiftY},{shiftX+1,shiftY}]

(* vrati pruseciky krusnize Circle[{u,v},r] s useckou x=[-1,1], y=0 *)
lineIntersect[u_,v_,radius_,shiftX_,shiftY_] :=
Module[{resList,root1,root2},
resList = {};
If[shiftY-radius≤v ≤shiftY+radius,
(* kdyz v je v intervalu shiftY-radius..shiftY+radius *)
(
(* spocitej pruseciky s primkou y=shiftY *)
root1 = u+Sqrt[(radius^2)-(((shiftY-v))^2)];
root2 = u-Sqrt[(radius^2)-(((shiftY-v))^2)];
(* zkontroluj jestli reseni lezi na usecce shiftX-1..shiftX+1 *)
If[-1+shiftX≤root1≤1+shiftX,AppendTo[resList,{root1,shiftY}]];
If[-1+shiftX≤root2≤1+shiftX,AppendTo[resList,{root2,shiftY}]];
)
];
Return[resList];
];

lineIntersect::usage = "lineIntersect[u_,v_,radius_,shiftX_,shiftY_] vrátí seznam bodů, ve kterých se protíná kružnice Circ[{u,v},radius] s úsečkou [{shiftX-1,shiftY},{shiftX+1,shiftY}]";

Průsečíky kružnice Circ[{u,v},radius] s obloukem Circ[{shiftX,shiftY},1,{0,Pi}]

(* vrati pruseciky krusnize Circle[{u,v},r] s pulkriznici Circle[{shiftX,shiftY},1,{0,Pi}] *)
arcIntersect[u_,v_,radius_,shiftX_,shiftY_] :=
Module[{curSol,allSol},
allSol = {};
(* spocitej pruseciky s kruznici Circ[{shiftX,shiftY},1] *)
curSol = DeleteDuplicates[NSolve[(((x-shiftX))^2)+(((y-shiftY))^2)== 1 && (((x-u))^2)+(((y-v))^2)== (radius^2),{x,y},Reals]];
(* pokud lezi v polorovine ≥shiftY, jsou to reseni pro pulkruznici *)
For[i=1,i≤Length[curSol],i++,If[curSol[[i,2,2]]≥ shiftY,AppendTo[allSol,{curSol[[i,1,2]],curSol[[i,2,2]]}]]];
Return[allSol];
];

arcIntersect::usage = "arcIntersect[u_,v_,radius_,shiftX_,shiftY_] vrátí průsečíky kružnice Circ[{u,v},radius] s půlkružnicí Circ[{shiftX,shiftY},1,{0,Pi}]";

Kruhová úseč

(* spocita kruhovou usec ze zadaneho uhlu a polomeru *)
usec[alpha_,r_] := (1/2)*(r^2) (alpha - Sin[alpha]);

usec::usage = "usec[alpha_,r_] vrátí plochu kruhové úseče kruhu se poloměrem r a úhlem mezi body an oblouku alpha (pouze konvexní úhel)";

(* spocita kruhovou usec ze zadanych body a polomeru -- zahrnuje i nekonvexni uhel, ale pouze specialne pro pripad pulkruhu a kruhu s mensim polomerem *)
usecAB[points_,u_,v_,r_] := 
Module[{alpha,result},
result = 0;
alpha =  Re[VectorAngle[SetPrecision[points[[1]]-{u,v},100],SetPrecision[ points[[2]] - {u,v},100]]];
result = usec[alpha,r];
If[(points[[1,2]]+ points[[2,2]])/2 ≥ v,result = Pi*(r^2)- result];
Return[result];
];

usecAB::usage = "usecAB[points_,u_,v_,r_] spočítá plochu úseče kruhu o poloměru r a středu {u,v} v závislosti na bodech na kruhovém oblouku, zohledňuje nekonvexní úhel, ale jen specificky (pokyd je v >= střed úsečky spojující body points -- konkávní)";

(* spocita kruhovou usec ze zadanych bodu a polomeru -- vzdy konvexni uhel mezi body points*)

usecABstatic[points_,u_,v_,r_] := 
Module[{alpha},
alpha =  Re[VectorAngle[SetPrecision[points[[1]]-{u,v},100],SetPrecision[ points[[2]] - {u,v},100]]];
Return[usec[alpha,r]];
];

usecABstatic::usage = "usecABstatic[points_,u_,v_,r_] spočítá plochu úseče kruhu o poloměru r a středu {u,v} v závislosti na bodech na kruhovém oblouku, vždy přacuje s konvexním úhlem mezi body points";

(* spocita prunik dvou kruznic ze zadanych bodu a polomeru -- prvni(pul)kriznice A *)
prunikAB[points_,u_,v_,r_,shiftX_,shiftY_] := 
Module[{result,segA,segB},
segA = usecABstatic[points,shiftX,shiftY,1];
segB = usecAB[points,u,v,r];
result = segA+segB;
Return[result];
];

prunikAB::usage = "prunikAB[points_,u_,v_,r_,shiftX_,shiftY_] Spočítá plochu průniku Disk[{u,v},r] a Disk[{shiftX,shiftY},1]. Předpokládá points jako dva body - průsečíky kružnic Circle[{u,v},r] a Circle[{shiftX,shiftY},1]";

Plocha průniku Disk[{u,v},r] a Disk[{shiftX,shiftY},1,{0,Pi}] = A

(* True pokud je stred {u,v} uvnitr mnoziny A *)
isCenterInside[u_,v_,shiftX_,shiftY_] :=
 If[shiftX-1≤ u≤shiftX+1 && Element[Sqrt[1-(((u-shiftX))^2)],Reals] &&  shiftY≤v≤shiftY+Sqrt[1-(((u-shiftX))^2)],True,False];

isCenterInside::usage = "isCenterInside[u_,v_,shiftX_,shiftY_] Vrátí true, pokud je bod {u,v} uvnitř Disk[{shiftX,shiftY},1,{0,Pi}]";

(* specialne pro pripad 1 pruseciku s Circ[{shiftX,shiftY},1,{0,Pi}] a úsečkou [{shiftX-1,shiftY},{shiftX+1,shiftY}] *)
intersectPart[arcInter_, lineInter_,u_,v_,r_,shiftX_,shiftY_] :=
Module[{arcX,arcY,lineX,segment, integral},
arcX = arcInter[[1,1]];
arcY = arcInter[[1,2]];
lineX = lineInter[[1,1]];

(* polovina úseče {arcX,-arcY} -- {arcX,arcY} kruhu a poloměrem 1 a středem {shiftX,shiftY} *)
segment=If[u≤  shiftX, 
If[lineX >  shiftX, (Pi/2) - (1/2)*usecABstatic[{{arcX,arcY},{arcX,shiftY-Abs[arcY-shiftY]}},shiftX,shiftY,1],(1/2)*usecABstatic[{{arcX,arcY},{arcX,shiftY-Abs[arcY-shiftY]}},shiftX,shiftY,1]],
If[lineX < shiftX,  (Pi/2) - (1/2)*usecABstatic[{{arcX,arcY},{arcX,shiftY-Abs[arcY-shiftY]}},shiftX,shiftY,1],(1/2)*usecABstatic[{{arcX,arcY},{arcX,shiftY-Abs[arcY-shiftY]}},shiftX,shiftY,1]]];

(* integrál přes zbylou oblast *)
If[Abs[shiftX-u] ≥ Abs[shiftX-arcX],
integral = If[u≤shiftX,
NIntegrate[Integrate[1,{x,arcX,Sqrt[(r^2)-(((y-v))^2)]+u}],{y,shiftY,arcY}],
NIntegrate[Integrate[1,{x,-Sqrt[(r^2)-(((y-v))^2)]+u,arcX}],{y,shiftY,arcY}]],
integral = If[u≤shiftX,
NIntegrate[Integrate[1,{y,shiftY,Sqrt[(r^2)-(((x-u))^2)]+v}],{x,arcX,lineX}],
NIntegrate[Integrate[1,{y,shiftY,Sqrt[(r^2)-(((x-u))^2)]+v}],{x,lineX,arcX}]]
];

(* oprava o úseč, ktera nebyla zapocitana kvuli x-ove mezi integralu *)
If[Abs[shiftX-u] <Abs[shiftX-arcX],
If[v>shiftY,
integral= integral +usecABstatic[{{lineX,shiftY},{lineX,shiftY+2Abs[shiftY-v]}},u,v,r];
]];

Return[segment+integral];
];

intersectPart::usage = "intersectPart[arcInter_, lineInter_,u_,v_,r_,shiftX_,shiftY_] Volána ve funkci intersect[] speciálně pro případ 1 pruseciku s Circ[{shiftX,shiftY},1,{0,Pi}] a úsečkou [{shiftX-1,shiftY},{shiftX+1,shiftY}]";

(* vraci obsah plochy pruniku A a Disk[{u,v},r] *)
intersect[u_,v_,r_,shiftX_,shiftY_] :=
Module[{arcInter, lineInter, arcNum, lineNum,result},
arcInter = arcIntersect[u,v,r,shiftX,shiftY];
lineInter = lineIntersect[u,v,r,shiftX,shiftY];
arcNum = Length[arcInter];
lineNum = Length[lineInter];

Switch[arcNum,
(* ZADNY PRUSECIK S OBLOUKEM *)
0, If[isCenterInside[u,v,shiftX,shiftY], 
(* stred jadra je UVNITR pulkruhu *)
Switch[lineNum,
0,result = Pi * (r^2),(* cele jadro je uvnitr *)
1,result = Pi * (r^2),(* cele jadro je uvnitr a dotyka se pulkruhu *)
2,result = Pi * (r^2) - usecAB[lineInter,u,v,r]], (* jadro je uvnitr castecne, usec je vne *)
(* stred jadra je VNE pulkruhu *)
Switch[lineNum,
0,result = 0,(* cele jadro je vne *)
1,result = 0,(* cele jadro je vne a dotyka se pulkruhu *)
2,result = usecABstatic[lineInter,u,v,r]] (* jadro je vne castecne, usec je uvnitr *)
],
(* 1 PRUSECIK S OBLOUKEM *)
1, If[isCenterInside[u,v,shiftX,shiftY], 
(* stred jadra je uvnitr pulkruhu *)
Switch[lineNum,
0,result = Pi * (r^2),(* cele jadro je uvnitr a dotyka se oblouku *)
1,result = intersectPart[arcInter, lineInter,u,v,r,shiftX,shiftY],(* PROBLEM *)
2,result = Pi * (r^2) - usecAB[lineInter,u,v,r]], (* jadro je uvnitr castecne, dotyka se oblouku, usec je vne *)
(* stred jadra je vne pulkruhu *)
Switch[lineNum,
0,result = 0,(* cele jadro je vne *)
1,result = intersectPart[arcInter, lineInter,u,v,r,shiftX,shiftY],(* PROBLEM *)
2,result = usecABstatic[lineInter,u,v,r]] (* jadro je vne castecne, usec je uvnitr *)
],
(* 2 pruseciky s obloukem *)
2, Switch[lineNum,
0,result = prunikAB[arcInter,u,v,r,shiftX,shiftY],(* cele jadro je uvnitr a dotyka se oblouku *)
1,result = prunikAB[arcInter,u,v,r,shiftX,shiftY],
2,If[isCenterInside[u,v,shiftX,shiftY], result = prunikAB[arcInter,u,v,r,shiftX,shiftY] - usecAB[lineInter,u,v,r], result = prunikAB[arcInter,u,v,r,shiftX,shiftY] - usecABstatic[lineInter,u,v,r]
]
];
];

(* Return[result]; *)
Return[1/(Pi*(r^2))*result];
];

intersect::usage = "intersect[u_,v_,r_,shiftX_,shiftY_] Vrací obsah plochy průniku A = Disk[{shiftX,shiftY},1,{0,Pi}] s Disk[{u,v},r]. Vypočteno numericky (přibližně)";

Alternativní výpočet plochy průniku A = Disk[{shiftX,shiftY},1,{0,Pi}] s Disk[{u,v},r] za pomoci RegionIntersection jádra Mathematica

(* urci obsah pruniku dvou kruznic -- vyuziva funkci jadra mathematicy -- RegionMeasure *)
regIntersect[x_,y_,diskR_,A_]:=
Module[{result,region},
region = RegionIntersection[A,Disk[{x,y},diskR]];
(*Print[DiscretizeRegion[region]];*)
If[RegionDimension[region] == 0, result = 0,result = RegionMeasure[region]];
(* Return[result];*)
Return[1/(Pi *(diskR^2))*result];
];

regIntersect::usage = "areaOfIntersectionOfAandDiskOfRadiusR[x_,y_,diskR_,A_] Určí obsah pruniku dvou kruznic -- vyuziva funkci jadra mathematicy -- RegionMeasure";

Hledání bodů na α-hladině

Hledání jednoho bodu v daném směru

(dirX,dirY):    představuje směr hledání,
height:          odpovídá výšce α  hledané α-hladiny,
epsilon > 0:  je maximální rozdíl mezi stupněm příslustnosti nalezeného bodu a height,
maxIter:        je maximální počet itrací

findLayerPoint[dirX_,dirY_,height_,epsilon_,maxIter_] :=
Module[{curX,curY,nextX,nextY, currentValue, nextValue, alpha,iterations,shiftX,shiftY},
If[height ≤ 0 || height > 1, Print["Hladina musí být mezi (0,1]"];Return[{}]];

{curX,curY} = {muStartX,muStartY}; (* pocatecni bod hledani (Subscript[x, 0]) *)
currentValue = fuzzyMemFunction[curX,curY]; (* pocatecni hodnota -- bude rovna 1 *)
iterations = 0; (* pocet iteraci *)
alpha = 1.; (* skalovaci faktor -- delka kroku *)


(* hledej dokud neprekrocis maximalni pocet iteraci nebo neni dosazeno pozadovane presnosti *)
While[Abs[currentValue-height]≥ epsilon && iterations < maxIter,
(* vypocti dalsi potencialni bod (Subscript[x, i+1]) *)
{nextX,nextY} = {curX,curY} + alpha * {dirX,dirY};


(* vypocti stupen prislusntosti v (Subscript[x, i+1]) *)
nextValue = fuzzyMemFunction[nextX,nextY];

(* pokud se priblizil cili a stupen prislusnosti je mensi -- presun se do nej a zvys alpha*)
If[nextValue ≥ height, 
(* kdyz ano *)
(
currentValue = nextValue;
{curX,curY} = {nextX,nextY};
alpha *= alphaSCALE;
),
(* kdyz ne, sniz alpha *)
alpha /= alphaSCALE;
];

iterations += 1;
];

(*Print["Počet iterací: ",iterations];*)
If[iterations==maxIter,Print["Pozor! Dosazeno maximalniho poctu iteraci ve smeru: ",{dirX,dirY}]];
Return[{curX,curY}];
];

findLayerPoint::usage = "findLayerPoint[dirX_,dirY_,height_,epsilon_,maxIter_]  najde bod, ve kterém je stupen prislusnosti fuzzy mnoziny (hodnota fuzzyMemFunction) priblizne roven pozadovane height s maximalni odchylkou epsilon, maxIter je maximalni pocet iteraci hledani, {dirX,dirY} smer hledani";

Hledání všech bodů na dané hladině

findAllLayerPoints[number_,height_,epsilon_,maxIter_] :=
Module[{directions,layerPoints,i,shiftX,shiftY},
If[height>1 || height ≤0,Print["Hladina musí být mezi (0,1]"];Return[{}]];

layerPoints = Parallelize[Map[findLayerPoint[#[[1]],#[[2]],height,epsilon,maxIter] &,N[generateCirclePoints[number], PRECISION] ]];
Return[layerPoints];
];
findAllLayerPoints::usage = "findAllLayerPoints[number_,height_,epsilon_,maxIter_] vrati seznam bodu fuzzy mnoziny na zadane hladine";

Hodnota opěrné funkce násobená směrem x/y

findSupport[{ex_,ey_},pts_] := Max[Map[({ex,ey} . #) &, pts]]
findSupport::usage="findSupport[{ex,ey},pts] najde maximum skalarniho soucinu smeru (ex,ey) a bodu z pts";

Hodnoty integrandu při výpočtu Steinerova bodu (hodnota opěrné funkce násobená souřadnicí směru):

supportTimesX[{ex_,ey_},pts_] := ex * findSupport[{ex,ey},pts];
supportTimesY[{ex_,ey_},pts_] := ey * findSupport[{ex,ey},pts];

Numerická integrace

Lichoběžníková metoda

Pro y-ovou souřadnici:

NtrapezoidY[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesY[#,N[pts[[;;;;Length[pts]/num ]],PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi; 
(* / (Pi * (1^2)) -- tj. obsah jednotkoveho kruhu -- konstanta pri vypoctu Steinerova bodu 2D *)

NtrapezoidY::usage="NtrapezoidY[num,pts] -- vygeneruje [num] uzlovych bodu na jednotkove kruznici a pocita y-ovou souradnici Steinerova bodu pomoci lichobeznikove metody; hodnoty operne funkce hledany pres pts (nikoliv celou A), delka pts musi byt mocninou dvou a minimalne maxPoints";

NtrapezoidYNoScale[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesY[#,N[pts,PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi;

Pro x-ovou souřadnici:

NtrapezoidX[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesX[#,N[pts[[;;;;Length[pts]/num ]],PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi; 
NtrapezoidX::usage="NtrapezoidY[num,pts] -- vygeneruje [num] uzlovych bodu na jednotkove kruznici a pocita x-ovou souradnici Steinerova bodu pomoci lichobeznikove metody; hodnoty operne funkce hledany pres pts (nikoliv celou A), delka pts musi byt mocninou dvou a minimalne maxPoints";

NtrapezoidXNoScale[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesX[#,N[pts,PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi;

Simpsonova metoda

NsimpsonY[num_,pts_] := 
(((2*Pi)/num)/3)* (1/ (Pi))*
Total[Table[If[Mod[i,2]==0,2,4],{i,0,num-1}] * 
ParallelMap[
(supportTimesY[ #,N[pts,PRECISION]]) &,N[generateCirclePoints[num],PRECISION]]];

Rombergova metoda

Při volání této metody je třeba, aby délka [pts] byla mocninou 2 a byla >= maxPoints, jinak je třeba volat NrombergYNoScale

NrombergY[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidY[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];
NrombergYNoScale[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidYNoScale[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

Stejně pro x-ovou souřadnici

NrombergX[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidX[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

NrombergXNoScale[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidXNoScale[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

Odhad chyby přes “horní” a “spodní” mnohoúhelník

Dělení bodů

V případě konvoluce s kruhem při použití množiny A vzniká “dole” úsečka (nalezená nepřesně). chceme nechat jen její krajní body a bod uprostřed

(* v pripade konvoluce A s kruhovym jaderem -- vznika dole usecka -- chceme odstranit prebytecne body *)
getLRMinOfList[points_] :=
Module[{i,list,minimum,left,right,flag},
list = points;
If[Length[list]<16,Print["Méně než 16. bodů."]];
minimum = list[[1,2]];
left = 0;
right = 0;
flag = False;


For[i=1,i<Length[list],i++,
If[(list[[i,2]]≤   minimum && Abs[list[[i,2]]-  minimum ] > (10^(-8))) ,
minimum = list[[i,2]];
left = 0;
right = 0;
flag = True;
];

If[ Abs[list[[i,2]]-  minimum ] ≤  (10^(-8)),flag = True];

If[flag == True,
If[list[[i,1]] < left, left = list[[i,1]]];
If[list[[i,1]] > right, right = list[[i,1]]];
flag = False;
]
];

Return[{{left,minimum},{right,minimum}}];
];
(* rozdeli [points] na dulezite body pro odhad a zbytek (na usecce) *)
getNewList[points_] :=
Module[{nums,numA,numB},
nums = getLRMinOfList[points];
numA = Position[points,nums[[1]]][[1,1]];
numB = Position[points,nums[[2]]][[1,1]];
Return[{
Join[points[[1;;numA]],{points[[(numA+numB)/2]]},points[[numB;;]]],
Join[points[[numA+1;;(numA+numB)/2-1]], points[[(numA+numB)/2+1;;numB-1]]]
}];
];

Opsaný mnohoúhelník L

<<ComputationalGeometry`

Funkce, která vrátí další/předchozí index (v kruhové smyčce)

nextPlus[number_,modulo_] := Mod[number,modulo]+1;
nextMinus[number_,modulo_] := If[number-1≤ 0,modulo,number-1];

Průsečík přímek daný body s indexy (i,i+1) a (i+2,i+3)

getIntersect[points_,i_]:=
Module[{aa,bb,cc,dd,AA,BB,CC,DD,mod},
mod=Length[points];
{aa,bb} = points[[i]];
{cc,dd} = points[[nextPlus[i,mod]]];
{AA,BB} = points[[ nextPlus[nextPlus[i,mod],mod] ]];
{CC,DD} = points[[ nextPlus[nextPlus[nextPlus[i,mod],mod],mod] ]];
 Return[{aa,bb}+(((-bb+BB) CC+AA (bb-DD)+aa (-BB+DD))/(-BB cc+bb CC-CC dd+AA (-bb+dd)+aa (BB-DD)+cc DD))*{aa-cc,bb-dd}];
];
getIntersect::usage="getIntersect[points,i] vrati prusecik primek dany body s indexy (i,i+1) a (i+2,i+3)";

getHull[allPoints_] :=
Module[{intersect,i,mod,hull,points,leftovers},
points = allPoints[[ConvexHull[allPoints]]];
(*{points,leftovers} = getNewList[allPoints];*)
(*points = allPoints;*)
mod = Length[points];
hull = {};

For[i=1,i≤ mod,i++,
intersect = getIntersect[points,i];
AppendTo[hull,intersect]
];

Return[hull];
];
getHull::usage="getHull[points] vrati vrcholy mnohouhelnika opsaneho mnohouhelniku urceneho vrcholy [points]";

Horní a spodní část

getUpperAndLower[points_] := Module[{hullOfPoints, mapping,iMaxRight,iMaxLeft, rightUpper, leftUpper,rightLower,leftLower, iRightUpper,iLeftUpper,iRightLower,iLeftLower, upperHull,lowerHull,upperPoint,lowerPoint,xCoordRight,xCoordLeft},
hullOfPoints = getHull[points];
mapping = Map[({1,0}. # ) &, points];
(* indexy bodu s maximalnim kolmym prumetem do (1,0) a (-1,0) z [points] *)
iMaxRight = Ordering[mapping,-1][[1]];
iMaxLeft = Ordering[mapping,1][[1]];

(* body a indexy pro opsany mnohouhelnik*)
xCoordRight = points[[iMaxRight]][[1]];
xCoordLeft = points[[iMaxLeft]][[1]];

(* iteruj pravy horni bod *)
While[True,
rightUpper =  getIntersect[points,nextMinus[iMaxRight,Length[points]]];
rightLower = getIntersect[points,nextMinus[nextMinus[iMaxRight,Length[points]],Length[points]]];
If[rightUpper[[1]] ≤ xCoordRight, Break[]];
iMaxRight = nextPlus[iMaxRight,Length[points]];
];

(* iteruj levy hodni bod*)
While[True,
leftUpper = getIntersect[points,nextMinus[nextMinus[iMaxLeft,Length[points]],Length[points]]];
leftLower =    getIntersect[points,nextMinus[iMaxLeft,Length[points]]];
If[leftUpper[[1]] ≥  xCoordLeft, Break[]];
iMaxLeft = nextMinus[iMaxLeft,Length[points]];
];

(*
rightUpper =  getIntersect[points,nextMinus[iMaxRight,Length[points]]];
leftUpper =    getIntersect[points,nextMinus[nextMinus[iMaxLeft,Length[points]],Length[points]]];
rightLower = getIntersect[points,nextMinus[nextMinus[iMaxRight,Length[points]],Length[points]]];
leftLower =    getIntersect[points,nextMinus[iMaxLeft,Length[points]]];
*)

iRightUpper = Position[hullOfPoints,rightUpper][[1,1]];
iLeftUpper =   Position[hullOfPoints,   leftUpper][[1,1]];
iRightLower = Position[hullOfPoints,rightLower][[1,1]];
iLeftLower =   Position[hullOfPoints,   leftLower][[1,1]];

upperHull = If[iRightUpper > iLeftUpper, 
Join[hullOfPoints[[iRightUpper;;Length[hullOfPoints]]], hullOfPoints[[1;;iLeftUpper]]],
Join[hullOfPoints[[iRightUpper;;iLeftUpper]]]
];
	lowerHull = If[iRightUpper > iLeftUpper, 
hullOfPoints[[iLeftUpper+1;;iRightUpper]][[1;;-2]],
Join[hullOfPoints[[iLeftUpper+1;;Length[hullOfPoints]]],hullOfPoints[[1;;iRightUpper]][[1;;-2]]]
];
	upperPoint =  If[iMaxRight > iMaxLeft, 
Join[points[[iMaxRight;;Length[points]]], points[[1;;iMaxLeft]]][[2;;-2]],
Join[points[[iMaxRight;;iMaxLeft]]][[2;;-2]]
];
	lowerPoint =  If[iMaxRight > iMaxLeft, 
points[[iMaxLeft+1;;iMaxRight]][[1;;-2]],
Join[points[[iMaxLeft+1;;Length[points]]],points[[1;;iMaxRight]][[1;;-2]]]
];

Return[{Join[upperHull,lowerPoint],Join[upperPoint,lowerHull]}] (* vrat sady bodu pro {horni,spodni} odhad *)
];

Vizualizace množin pro horní a spodní odhad

visualiseHull[points_] := Module[{hullOfPoints,iMaxRight,iMaxLeft,upperSet,lowerSet,mapping},
hullOfPoints = getHull[points];
{upperSet,lowerSet} = getUpperAndLower[points];

mapping = Map[({1,0}. # ) &, points];
iMaxRight = Ordering[mapping,-1][[1]];
iMaxLeft = Ordering[mapping,1][[1]];

Print[
Graphics[{
Black,
Point[points],
Point[hullOfPoints],
Line[Append[points,points[[1]]]],
Line[Append[hullOfPoints,hullOfPoints[[1]]]] 
}]
];

Print[
Graphics[{
PointSize[Small],
Black,
Point[points],
Point[hullOfPoints],
Line[Append[points,points[[1]]]],
(*Line[Append[hullOfPoints,hullOfPoints[[1]]]] ,*)
Text[Subscript[k,"+"],points[[iMaxRight]] - {0.06,0}],
Text[Subscript[k,"-"],points[[iMaxLeft]]+{0.06,0}],
Red, 
Line[Append[upperSet,upperSet[[1]]]],
Blue,
Line[Append[lowerSet,lowerSet[[1]]]]
}]
];
]

Hledání hladin, Steinerových bodů, vizualizace ...

Hledání hladin

findAndSaveLayers[heights_,heightsLabels_,numberOfPoints_, epsilon_,maxIter_]:=
Module[{index,data},
For[index=1,index≤Length[heights],index++,
Print["Current index: ",index, ", height ", heights[[index]], ", filename: ",heightsLabels[[index]]];
data = findAllLayerPoints[numberOfPoints,heights[[index]],epsilon,maxIter];
WriteString["stdout","\n"];
Export[heightsLabels[[index]],data];
];
]
findAndSaveLayers::usage = "findAndSaveLayers[heights,heightsLabels,numberOfPoints, epsilon,maxIter] najde [numberOfPoints] bodu na kazde z [heights] hladin s presnosti [epsilon], ulozi je do soboru [heightsLabels], pri hledani kazdeho bodu nesmi prekrocit [maxIter], jinak vypise varovani";

Počítání Steinera

computeAndSaveSteiner[heightsLabels_,steinerLabels_] :=
Module[{layerData,layerSteiner,romTable,steinerUpperEst,steinerLowerEst,estTable,upperSet,lowerSet,i,j,tALower,tAUpper},
layerData = Array[0 &,Length[heightsLabels]];
layerSteiner = {};steinerUpperEst = {}; steinerLowerEst= {};
Print["Load data and compute Steiner -- start"];

Print["  -- Importing data"];
(* nacti data *)
For[i=1,i≤Length[heightsLabels],i++, layerData[[i]] = Import[heightsLabels[[i]]]];

Print["  -- Computing Steiner points"];
(* spocitej steinerovy body *)
For[i=1,i≤Length[layerData],i++, 
romTable = integrationMETHOD[maxPoints,layerData[[i]] ];
AppendTo[layerSteiner,Diagonal[romTable] ];

(* spocitej odhady  *)
tALower = {}; tAUpper={};
For[j=minPoints,j≤maxPoints,j*=2,
{upperSet,lowerSet} = getUpperAndLower[layerData[[i]][[;;;; Length[layerData[[i]]] /j ]] ];
estTable = integrationMETHODNS[maxPoints,upperSet];
AppendTo[tAUpper,estTable[[Length[estTable],Length[estTable]]] ];

estTable = integrationMETHODNS[maxPoints,lowerSet];
AppendTo[tALower,estTable[[Length[estTable],Length[estTable]]] ];
];

AppendTo[steinerUpperEst,tAUpper];
AppendTo[steinerLowerEst,tALower];
];

Print["  -- Exporting results"];
(* uloz Steinerovy body *)
For[i=1,i≤Length[layerSteiner],i++, Export[steinerLabels[[i]],{layerSteiner[[i]] , steinerUpperEst[[i]],steinerLowerEst[[i]]}]
];
];
computeAndSaveSteiner::usage = "computeAndSaveSteiner[heightsLabels,steinerLabels] spocita steinerovy body z bodu na hladinach, ktere nacte ze souboru [heightsLabels] a ulozi je do souboru [steinerLabels]";

Vizualizace

printVisuals[heights_,heightsLabels_,steinerLabels_] :=
Module[{layerData,layerSteiner,i,layerData3D,plot3Dlist,plot3Dorig, layerSteinerAll,contours,xMin,xMax,yMin,yMax, layersGraphis,layerSteinerUpper,layerSteinerLower,importBuf,upperEsts,lowerEsts,estGraphic},
layerData = Array[0 &,Length[heightsLabels]];
layerSteiner = {};layerSteinerUpper={}; layerSteinerLower={};upperEsts={}; lowerEsts={};estGraphic={};
Print["Load data and compute Steiner -- start"];

Print["  -- Importing data"];
(* load layer data *)
For[i=1,i≤Length[layerData],i++, layerData[[i]] = Import[heightsLabels[[i]]]];
For[i=1,i≤Length[steinerLabels],i++,
importBuf = Import[steinerLabels[[i]] ];
AppendTo[layerSteiner,importBuf[[1]]];
AppendTo[layerSteinerUpper,importBuf[[2]]];
AppendTo[layerSteinerLower,importBuf[[3]]];

AppendTo[upperEsts,Flatten/@Transpose[{Array[#&,Length[layerSteinerUpper[[i]] ]],layerSteinerUpper[[i]]}] ];
AppendTo[lowerEsts,Flatten/@Transpose[{Array[#&,Length[layerSteinerLower[[i]] ]],layerSteinerLower[[i]]}] ];

];

	(* plot range *)
xMin = MinimalBy[layerData[[1]],First][[1,1]];
       xMax = MaximalBy[layerData[[1]],First][[1,1]];
       yMin = MinimalBy[layerData[[1]],Last][[1,2]];
       yMax = MaximalBy[layerData[[1]],Last][[1,2]];

Print["  -- Visualisation"];

Print["     :plot all layers"];
Print[ListPlot[layerData,AspectRatio->Automatic,PlotStyle->Map[ColorData["Rainbow"][#]&,Prepend[Array[1/(Length[layerData]-1)*#&,Length[layerData]-1],0]]]];

Print["     :plot layers in 3D"];
layerData3D = Array[0 &,Length[layerData]];
	For[i=1,i≤Length[layerData],i++,
layerData3D[[i]] = MapThread[Append,{layerData[[i]],Array[heights[[i]] &, Length[layerData[[i]] ] ]}]
];

plot3Dlist =ListPointPlot3D[Flatten[layerData3D,1],AxesLabel->{"x",Framed["y",FrameStyle->None,FrameMargins->10],Framed["z",FrameStyle->None,FrameMargins->10]}, ColorFunction->"Rainbow",PlotRange->All,AspectRatio->Automatic];
Print[plot3Dlist];
(*
plot3Dorig =Plot3D[regIntersect[x,y,R,setA],{x,xMin,xMax},{y,yMin,yMax},PlotRange->All,AxesLabel->{"x",Framed["y",FrameStyle->None,FrameMargins->10],Framed["z",FrameStyle->None,FrameMargins->10]},AspectRatio->Automatic];
Print[Show[{plot3Dlist,plot3Dorig}]];

Print[plot3Dorig];


Print["     :contours"];
contours = ContourPlot[regIntersect[x,y,R,setA],{x,xMin,xMax},{y,yMin,yMax},Contours->heights,AspectRatio->Automatic,PlotLegends->Automatic];
Print[contours];
Print[Show[{contours, ListPlot[layerData,AspectRatio->Automatic]}]];
*)

Print["     :visualize Steiner points"];
layerSteinerAll = {};
For[i=1,i≤ Length[layerSteiner],i++,AppendTo[layerSteinerAll, layerSteiner[[i,-1]]];];

Print[layerSteinerAll];

Print[ListPlot[layerSteinerAll,Frame->True,GridLines->Automatic,FrameLabel->{{"(Subscript[s, y])",None},{None,"α"}}, FrameTicks->{{Automatic,None},{None,Transpose[{Array[# &, Length[heights]],Map[N[#] &,heights]}] }} ]];

Print["     :table of Steiners"];
Print[TableForm[Transpose[layerSteiner],TableHeadings->{
Array[StringJoin["Bodů: ",ToString[2^(#+Log[2,minPoints]-1)]]&,Log[2,maxPoints]-Log[2,minPoints]+1],Map[StringJoin["α = ", ToString[N[#]]] &,heights]
}]];

Print["     :error estimates"];

For[i=1,i≤Length[layerSteiner],i++,
AppendTo[estGraphic,
Show[{
ListPlot[layerSteiner[[i]],PlotStyle->{Black,PointSize[Medium]},Ticks->{Flatten/@Transpose[{Array[#&,Length[Array[ToString[2^(#+Log[2,minPoints]-1)]&,Log[2,maxPoints]-Log[2,minPoints]+1]]],Array[ToString[2^(#+Log[2,minPoints]-1)]&,Log[2,maxPoints]-Log[2,minPoints]+1]}] ,Automatic}],
Graphics[{
PointSize[Medium],
Red,
Point[upperEsts[[i]]],
Blue,
Point[lowerEsts[[i]]]
}],
ListPlot[layerSteiner[[i]],PlotStyle->{Black,PointSize[Medium]}]
},PlotRange->All,PlotLabel->ToString[heightsLabels[[i]]],AxesLabel->{N,(Subscript[s, y])}]
];
];

Print[estGraphic];

Print["     :additional graphics"];

layersGraphis = {};
For[i=1,i≤  Length[layerData],i++,
AppendTo[layersGraphis,Graphics[{
Thickness[0.01],
Line[Append[layerData[[i]],layerData[[i,1]]]],
Gray,
Point[{{0,0}}],
PointSize[Large],
Red,
Point[{0,layerSteinerAll[[i]]}],
Black
(*Inset[ StringJoin["Hladina ", ToString[N[layerHeights[[i]]]] ] ,{0,0.05}]*)
},AspectRatio->Automatic,Axes->True, PlotLabel->StringJoin["Hladina ", ToString[N[heights[[i]]]], ", sy = ", ToString[layerSteinerAll[[i]]]],PlotRange->{{xMin-0.1,xMax+0.1},{yMin-0.1,yMax+0.1}}]
];
];

Print[layersGraphis]
];

printVisuals::usage="printVisuals[heights,heightsLabels,steinerLabels] provede kompletni vizualizaci - body na hladinach a vypoctene steinerovy body nacte se souboru [heightsLabels] resp. [steinerLabels]";

Funkce, která najde body, spočítá steinerovy body, exportuje a vizualizuje

doItAll[heights_,heightsLabels_,steinerLabels_,numberOfPoints_, epsilon_,maxIter_] :=
Module[{},
findAndSaveLayers[heights,heightsLabels,numberOfPoints, epsilon,maxIter];
computeAndSaveSteiner[heightsLabels,steinerLabels];
printVisuals[heights,heightsLabels,steinerLabels];
];
doItAll::usage="doItAll[heights,heightsLabels,steinerLabels,numberOfPoints, epsilon,maxIter] vola funkce findAndSaveLayers, computeAndSaveSteiner, printVisuals s prislusnymi parametry";