Scaling at depth 1	Scaling at depth 2

2. По-важни детайли.
1. За да бъде една окръжност с достатъчно голям радиус поместена в опашката като "хубава", тя, освен всичко друго трябва да бъде максимална в смисъл, че не може да се увеличи нейния радиус с малки деформации, такива че по време на тази деформация нито един връх на многоъгълника не влиза вътре в окръжността. Дефиниция на това, какво точно значи "малка деформация" е "дълга и широка", но интуитивно е ясно, че такива неразширими окръжности са в точност тези, които са породени от точки, образуващи остроъгълен триъгълник. Трябва да се включи и граничен случай на правоъгълните триъгълници (иначе ще бъдат игнорирани правоъгълниците!). А тъй като всички сравнения от типа X <= Y трябва да се провеждат "с точност до epsilon" (за разлика от изчисления с плаваща запетая, провеждани с максимална възможна тожност), някое число от "почти правоъгълни" триъгълници с тъп ъгъл попадат в съответните приоритетни опашки. Което е довело до един мръсен бъг, присъстващ в тестваната версия: "изчезването" на някои от такива окръжности от приоритетна опашка при двойно намаляване на параметра scale. Безуспешното търсене на причините за този бъг (който ме е лишил от 2 точки) ми е отнело почти една седмица по време на писането на кода, а едноредов bug fix го открих само преди няколко дена: ред 413 от моя сорс "fire.cpp" трябва да се замени с
return MaxCircle(depth-1,c0,scale/2);//"No circle" bug fix (30.3.02)
2. Друг важен момент: коректна работа с некоректните "плаващи" числа. Имаше едно явление, което ме тормозеше много: при изчисляването (по известни формули от училището) на цънтъра на окръжността породена от един триъгълник с много остър ъгъл грешката е била много голяма -- някъде в 4-ти или 5-ти знак след запетайката (дори при изчисления с long double!), което при по-малки epsilon пак довеждаше до неприятни явления с изчезването на окръжности при увеличаване на дълбочината на рекурсия. Ще напомня, че epsilon, параметър на точността на сравнения; при неговото увеличаване до размера на грешките при изчисляване на центъра на окръжността, броя на окръжности в опашките се увеличаваше многократно, драстично намалявайки скоростта на алгоритма.
Тъй че за по-точното намиране на центровете аз използвах след "училищните" изчисления и една quicky (в смисъл, написана бързо на ръка :) версия на метода на последователни приближения. Разбира се по-добре би било вместо това просто да се търси пресечката на симетрала на най-късата страна на триъгълника с симетрала на най-дългата му страна, защото големите грешки при малки ъгли възникват тъкмо заради това, че се търси пресечката на две почти успоредни линии. Обаче заради терапаната (писах този коригиращ код в последния ден вечерта, няколко часа преди изтичане на срока) написах това, което написах.
3. Изпъкналостта на многоъгълника играе в този алгоритм много важна роля. Тя позволява да прецени за една перспективна окръжност C с радиус radius (и построена от "разредени" точки), долната граница inf(radius) на радиусa на максималната "хубава" окръжност която може да се построи от всички (т.е. неразредени) точки на многоъгълника съдържащи се вътре в C. Така че, ако inf(radius) е по-голямо (с точност epsilon) от радиус на окръжността в главата на опашката, опашката се изпразва, след което в нея се помества окръжност C. Съответния параметр delta за намиране на тази долна граница (inf(radius)=radius-delta) аз го изчислих по една приближена формула delta=scale*scale/(radius*2) . Тази формула е добра, когато ъгъла delta=scale/radius е малък (какъвто е масовия случай), обаче понякога той може да се окаже голям, както аз открих преди няколко дена! Подозирам, че това е причина поне за някои от оставащи все още бъгове. А точната формула (с sin и cos) беше нетрудно да се напише, но ми беше мързел (и все още е :).
5. Коментари към "авторското" решение.
1. Авторското решение с hill climbing, т.е. решението описано във авторския анализ на задачата, възможно би могло да се приложи успешно за намиране на максимума (минимума) на непрекъсната функция f ако тази функция би била вдлъбната (изпъкнала), или, по-общо, би имала само един локален максимум(минимум), който едновременно е и глобален. Тъкмо такъв бе случая с 6-та задача от миналогодишния конкурс [1], където аз използвах една по-канонична разновидност на този метод (в смисел че е описана във всички учебници :) известна под кодово име subgradient method. При този метод: началната точка P се избира горе-долу произволно, но желателно да е в околността на центъра на тежестта (както и в авторското решение); вместо М точки се избират М единични вектори (направления) n1,...nM. След това се избира направление ni, такова че разликата d=f(P+epsilon*ni)-f(P) е максимална (минимална). Тук epsilon е исканата точност. Ако модула на това максимално d е по-малък от epsilon, значи P е оптимална с зададена точност. В противен случай намираме точката P' на права, минаваща през P в направление ni на най-бързото качване (спускане), в която функцията f, а по-точно нейното ограничение върху правата е оптимално. Полагаме P=P' и продължаваме цикла. Максимума върху правата се намира както е описано в авторското решение, с тези разлики че а) възможните направления в едномерния случай са само две и б)вместо с "голяма", но невъзрастваща (както е в авторското решение) стъпка S, се започва с "безкрайно малка" стъпка S=epsilon, като стъпката всеки следващ път се удвоява (или, по-добре се умножава със "златното сечение" 1.618), докато е възможно. Тогава интервал между последната и предпоследната точки на правата съдържи оптималната за тази права точка, която след това стесняваме този интервал с дихотомии: разделяме го на два равни подинтервала (или, по-добре, в отношение 1:1.618 -- пак "златното сечение"), определяме, в кой от тях се намира оптималната точка, и т.н., докато модул на разликата на целевата функция |f(P0)-f(P1)| в началото и в края на интервала (не дължината на интервала!) ще стане по-малък от epsilon.
Тези малки разлики довеждат до това че subgradient method е със сложност O(log(D/epsilon)) и в най-лошия случай, докато hill climbing гарантира само O(D/epsilon) (D е разстоянието между началната и оптималната точки). Да речем, ако областта, в която се търси оптимум е правоъгълник 1М*1, тогава не може да се избере начално S по-голямо от 1. И ако 1М-2 точки се намират на левия край на този дълъг правоъгълник близо една до друга и само две точки -- в ъглите на десния край, то центъра на тежестта се намира около ляв край, оптимална точка -- по средата, (като на рисунка: ) и алгоритм ще трае 1М*0.5М=5*10^11*средно време на елементарна операция T.e. 500 сек. на 1G процесор, ако средна елементарна операция се извершва за един такт -- което е далеч от реалността :-).

2. Обаче всичко казано горе важи само за хубавите функции, с единствен локален и глобален максимум. Колко такива локални максимуми може да има в нашата задача? Отговор: число от порядъка O(N): лесно се вижда, че в термините на диаграмите на Вороной, всеки локален минимум се намира във един от върховете на диаграмата на Вороной. Обратното не е вярно: в тъпоъгълен триъгълник единствен връх на диаграмата на Вороной (пресечна точка на тритте симетрали) не е точка на локален максимум, т.е. за такъв триъгълник локален максимум изобщо няма. Тъкмо това бе и причината, заради която аз изключих в алгоритма ми такива триъгълници при намирането на "вътрешна" оптимална точка. А сега примера, илюстриращ случая, когато локалните максимуми са много: (червените точки на картинката). От тази картинка се вижда също така, че както "авторския", така и "каноничен" вариант на метода на най-бързото качване яко ще се издъни на този простичек пример: точката на локалния максимум в която алгоритм ще се приземи силно зависи от точката на старта му. Ако се стартира от "цънтъра на тежестта" (синята точка), ще бъде върната като оптимална една от дветте съседни червени точки, докато същинския оптимум очевидно е най-лявата такава. Вижда се също така, че условието за единственост на оптималната точка по никакъв начин не помага в намирането на точка от Z-околността на оптималната, тъй че спокойно може да бъде махнато от формулировката на задачата. Нормален безбъгов алгоритм чисто и просто би трябвало да връща произволна точка, оптимизираща целевата функция със зададена точност epsilon [7].

3. Разбира се, при многократно прилагане на subgradient или hill climbing алгоритм, със случайно (и достатъчно голямо) множество от начални точки вместо една единствена, шансовете за намиране на коректен локален максимум (т.е. максимален измежду всичките локалните максимуми) се увеличават. Това става обаче за сметка на скоростта на алгоритма. Т.е. Монте-Карло версия на алгоритма макар и по-коректна, нада ли ще се окаже достатъчно ефективна.
6. Коментари към коментари.
[1] Aз приложих модификация на част от този алгоритм за решаването на 6-та задача на миналогодишния конкурс (Circles) която е дуална на разглежданата тук: там целевата функция f беше дефинирана като f(P)=max{distance(P,Pi)} и се търсеше нейния минимум, което драстично упростяваше нещата, заради това че f е изпъкнала в този случай, и за нейно изчисление е било достатъчно да се използват точки Pi само от изпъкналата обвивка на множеството от всички точки. А изпъкналата обвивка беше до стотици пъти по-малка от множеството на всички точки, заради спецификата на задачата. Е, загубих доста точки тогава заради използването на скапания scanf() (мизерен максимален буфер на който -- 32К -- е останал, MS знае защо, от епохата на DOS) вместо по-адекватен за големите файлове fread(). Всъщност аз не използвах оригиналния сорс на Qhull, а написах своя код (само за 2D) въз основа на авторското описание на Qhull "The Quickhull Algorithm for Convex Hulls". Кода ми беше 2 пъти по-бърз от Qhull, и аз бях много горд от това до момента, когато видях резултати и кода на Камен Добрев за изчисляване на изпъкнала обвивка: той беше 15 до 30 пъти по-бърз от моя! Този свръх бърз код исползваше до дупка специфични особености на задачата (целочислеността на координатите на точките и относително малките размери на областта, в която те се намират: квадрат 100000 на 100000), но не е изключено, че може да бъде обобщен и за по-общ случай. Както и да е, К.Добрев загуби тогава точки, защото негов brute force код за намиране на оптималната точка беше много бавен в случая, когато изпъкналата обвивка бе голяма (няколко хиляди точки). Впрочем имаше и един бъг: масив за точки от изпъкналата обвивка се е оказал малък за един от тестовите файлове.
Друг интересен момент в тази задача беше, че за нея съществува и друго, "native", решение, което не използва изпъкнали обвивки изобщо, намирайки директно вместо това нещо, което може да се нарече "кръгова обвивка". Този алгоритъм, намерен от Мартин Вълканов и Петър Петров, е много прост и бърз (те заеха първите 2 места тогава).

[2] Оригинална програма qvoronoi.exe е компилиран с много стара версия на Borland C++ от 1995 г. Нейната рекомпилация с VC++6 би трябвало да ускори програмата със около 20%.

[3] "Във много знание -- много печал", както е отбелязал Еклесиаст на времето. Разбира се обратното твърдение също е вярно, въпреки законите на логиката. Да речем, този, който знаеше ключовата дума "perfect match" (без дори да е чувал преди това за Blossom IV или V ;-) е имал шанс да спечели доста екстра точки от 3-та задача.

[4] Тъй че твърдението от "авторския" коментар че външната оптимална точка "...ще е или някой от върховете на многоъгълника или пресечна точка на контура със симетрала на страна на многоъгълника" не е вярно, според мен, в общия случай. А твърдението за O(N)-сложност на задачата при исползване на алгоритма на R.Klein и A.Lingas за исчисляване на ограничени диаграми на Вороной -- ако и е вярно, то поне не е толкова очевидно -- поне за мен.
Коментар на журито: Цитатът по-горе е грешен. Нашият коментар казва "В случай, че търсената точка е извън многоъгълника, не е трудно да се докаже, че това ще бъде точка, лежаща на страна на правоъгълника и това ще е или някои от върховете му или пресечна точка на контура със симетрала на страна на многоъгълника."

[5] Пример на един такъв лош случай ("змия"):

[6] Току-що открих и тест-убиец на scaling алгоритм ("двe змии" при N==34):