cd.lua 39 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532
  1. -- This code is derived from the SOM benchmarks, see AUTHORS.md file.
  2. --
  3. -- Copyright (c) 2016 Francois Perrad <[email protected]>
  4. --
  5. -- Permission is hereby granted, free of charge, to any person obtaining a copy
  6. -- of this software and associated documentation files (the 'Software'), to deal
  7. -- in the Software without restriction, including without limitation the rights
  8. -- to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  9. -- copies of the Software, and to permit persons to whom the Software is
  10. -- furnished to do so, subject to the following conditions:
  11. --
  12. -- The above copyright notice and this permission notice shall be included in
  13. -- all copies or substantial portions of the Software.
  14. --
  15. -- THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  16. -- IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17. -- FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  18. -- AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  19. -- LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  20. -- OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  21. -- THE SOFTWARE.
  22. --[[
  23. The module 'bit' is available with:
  24. * LuaJIT
  25. * LuaBitOp extension which is available for:
  26. * Lua 5.1
  27. * Lua 5.2
  28. The module 'bit32' is available with:
  29. * Lua 5.2
  30. * Lua 5.3 when compiled with LUA_COMPAT_5_2
  31. The bitwise operators are added to Lua 5.3 as new lexemes (there causes
  32. lexical error in older version)
  33. --]]
  34. local band, bxor, rshift
  35. if _VERSION < 'Lua 5.3' then
  36. local bit = bit32 or require'bit'
  37. band = bit.band
  38. bxor = bit.bxor
  39. rshift = bit.rshift
  40. else
  41. band = assert(load'--[[band]] return function (a, b) return a & b end')()
  42. bxor = assert(load'--[[bxor]] return function (a, b) return a ~ b end')()
  43. rshift = assert(load'--[[rshift]] return function (a, b) return a >> b end')()
  44. end
  45. local alloc_array
  46. local ok, table_new = pcall(require, 'table.new') -- LuaJIT 2.1 extension
  47. if ok then
  48. alloc_array = function (n)
  49. local t = table_new(n, 1)
  50. t.n = n
  51. return t
  52. end
  53. else
  54. alloc_array = function (n)
  55. local t = {}
  56. t.n = n
  57. return t
  58. end
  59. end
  60. local Vector = {_CLASS = 'Vector'} do
  61. local floor = math.floor
  62. function Vector.new (size)
  63. local obj = {
  64. storage = alloc_array(size or 50),
  65. first_idx = 1,
  66. last_idx = 1,
  67. }
  68. return setmetatable(obj, {__index = Vector})
  69. end
  70. function Vector.with (elem)
  71. local v = Vector.new(1)
  72. v:append(elem)
  73. return v
  74. end
  75. function Vector:at (idx)
  76. if idx > self.storage.n then
  77. return nil
  78. end
  79. return self.storage[idx]
  80. end
  81. function Vector:at_put (idx, val)
  82. if idx > self.storage.n then
  83. local new_n = self.storage.n
  84. while idx > new_n do
  85. new_n = new_n * 2
  86. end
  87. local new_storage = alloc_array(new_n)
  88. for i = 1, self.storage.n do
  89. new_storage[i] = self.storage[i]
  90. end
  91. self.storage = new_storage
  92. end
  93. self.storage[idx] = val
  94. if self.last_idx < idx + 1 then
  95. self.last_idx = idx + 1
  96. end
  97. end
  98. function Vector:append (elem)
  99. if self.last_idx > self.storage.n then
  100. -- Need to expand capacity first
  101. local new_storage = alloc_array(2 * self.storage.n)
  102. for i = 1, self.storage.n do
  103. new_storage[i] = self.storage[i]
  104. end
  105. self.storage = new_storage
  106. end
  107. self.storage[self.last_idx] = elem
  108. self.last_idx = self.last_idx + 1
  109. end
  110. function Vector:is_empty ()
  111. return self.last_idx == self.first_idx
  112. end
  113. function Vector:each (fn)
  114. for i = self.first_idx, self.last_idx - 1 do
  115. fn(self.storage[i])
  116. end
  117. end
  118. function Vector:has_some (fn)
  119. for i = self.first_idx, self.last_idx - 1 do
  120. if fn(self.storage[i]) then
  121. return true
  122. end
  123. end
  124. return false
  125. end
  126. function Vector:get_one (fn)
  127. for i = self.first_idx, self.last_idx - 1 do
  128. local e = self.storage[i]
  129. if fn(e) then
  130. return e
  131. end
  132. end
  133. return nil
  134. end
  135. function Vector:remove_first ()
  136. if self:is_empty() then
  137. return nil
  138. end
  139. self.first_idx = self.first_idx + 1
  140. return self.storage[self.first_idx - 1]
  141. end
  142. function Vector:remove (obj)
  143. local new_array = alloc_array(self:capacity())
  144. local new_last = 1
  145. local found = false
  146. self:each(function (it)
  147. if it == obj then
  148. found = true
  149. else
  150. new_array[new_last] = it
  151. new_last = new_last + 1
  152. end
  153. end)
  154. self.storage = new_array
  155. self.last_idx = new_last
  156. self.first_idx = 1
  157. return found
  158. end
  159. function Vector:remove_all ()
  160. self.first_idx = 1
  161. self.last_idx = 1
  162. self.storage = alloc_array(self:capacity())
  163. end
  164. function Vector:size ()
  165. return self.last_idx - self.first_idx
  166. end
  167. function Vector:capacity ()
  168. return self.storage.n
  169. end
  170. function Vector:sort (fn)
  171. -- Make the argument, block, be the criterion for ordering elements of
  172. -- the receiver.
  173. -- Sort blocks with side effects may not work right.
  174. if self:size() > 0 then
  175. self:sort_range(self.first_idx, self.last_idx - 1, fn)
  176. end
  177. end
  178. function Vector:sort_range (i, j, fn)
  179. assert(fn)
  180. -- The prefix d means the data at that index.
  181. local n = j + 1 - i
  182. if n <= 1 then
  183. -- Nothing to sort
  184. return
  185. end
  186. local storage = self.storage
  187. -- Sort di, dj
  188. local di = storage[i]
  189. local dj = storage[j]
  190. -- i.e., should di precede dj?
  191. if not fn(di, dj) then
  192. local tmp = storage[i]
  193. storage[i] = storage[j]
  194. storage[j] = tmp
  195. local tt = di
  196. di = dj
  197. dj = tt
  198. end
  199. -- NOTE: For DeltaBlue, this is never reached.
  200. if n > 2 then -- More than two elements.
  201. local ij = floor((i + j) / 2) -- ij is the midpoint of i and j.
  202. local dij = storage[ij] -- Sort di,dij,dj. Make dij be their median.
  203. if fn(di, dij) then -- i.e. should di precede dij?
  204. if not fn(dij, dj) then -- i.e., should dij precede dj?
  205. local tmp = storage[j]
  206. storage[j] = storage[ij]
  207. storage[ij] = tmp
  208. dij = dj
  209. end
  210. else -- i.e. di should come after dij
  211. local tmp = storage[i]
  212. storage[i] = storage[ij]
  213. storage[ij] = tmp
  214. dij = di
  215. end
  216. if n > 3 then -- More than three elements.
  217. -- Find k>i and l<j such that dk,dij,dl are in reverse order.
  218. -- Swap k and l. Repeat this procedure until k and l pass each other.
  219. local k = i
  220. local l = j - 1
  221. while true do
  222. -- i.e. while dl succeeds dij
  223. while k <= l and fn(dij, storage[l]) do
  224. l = l - 1
  225. end
  226. k = k + 1
  227. -- i.e. while dij succeeds dk
  228. while k <= l and fn(storage[k], dij) do
  229. k = k + 1
  230. end
  231. if k > l then
  232. break
  233. end
  234. local tmp = storage[k]
  235. storage[k] = storage[l]
  236. storage[l] = tmp
  237. end
  238. -- Now l < k (either 1 or 2 less), and di through dl are all
  239. -- less than or equal to dk through dj. Sort those two segments.
  240. self:sort_range(i, l, fn)
  241. self:sort_range(k, j, fn)
  242. end
  243. end
  244. end
  245. end -- class Vector
  246. local Set = {_CLASS = 'Set'} do
  247. local INITIAL_SIZE = 10
  248. function Set.new (size)
  249. local obj = {
  250. items = Vector.new(size or INITIAL_SIZE)
  251. }
  252. return setmetatable(obj, {__index = Set})
  253. end
  254. function Set:size ()
  255. return self.items:size()
  256. end
  257. function Set:each (fn)
  258. self.items:each(fn)
  259. end
  260. function Set:has_some (fn)
  261. return self.items:has_some(fn)
  262. end
  263. function Set:get_one (fn)
  264. return self.items:get_one(fn)
  265. end
  266. function Set:add (obj)
  267. if not self:contains(obj) then
  268. self.items:append(obj)
  269. end
  270. end
  271. function Set:remove_all ()
  272. self.items:remove_all()
  273. end
  274. function Set:collect (fn)
  275. local coll = Vector.new()
  276. self:each(function (it)
  277. coll:append(fn(it))
  278. end)
  279. return coll
  280. end
  281. function Set:contains (obj)
  282. return self:has_some(function (it) return it == obj end)
  283. end
  284. end -- class Set
  285. local IdentitySet = {_CLASS = 'IdentitySet'} do
  286. setmetatable(IdentitySet, {__index = Set})
  287. function IdentitySet.new (size)
  288. local obj = Set.new(size)
  289. return setmetatable(obj, {__index = IdentitySet})
  290. end
  291. function IdentitySet:contains (obj)
  292. return self:has_some(function (it) return it == obj end)
  293. end
  294. end -- class IdentitySet
  295. local Entry = {_CLASS = 'Entry'} do
  296. function Entry.new (hash, key, value, next)
  297. local obj = {
  298. hash = hash,
  299. key = key,
  300. value = value,
  301. next = next,
  302. }
  303. return setmetatable(obj, {__index = Entry})
  304. end
  305. function Entry:match (hash, key)
  306. return self.hash == hash and self.key == key
  307. end
  308. end -- class Entry
  309. local Dictionary = {_CLASS = 'Dictionary'} do
  310. local INITIAL_CAPACITY = 16
  311. function Dictionary.new (size)
  312. local obj = {
  313. buckets = alloc_array(size or INITIAL_CAPACITY),
  314. size = 0,
  315. }
  316. return setmetatable(obj, {__index = Dictionary})
  317. end
  318. function Dictionary:hash (key)
  319. if not key then
  320. return 0
  321. end
  322. local hash = key:custom_hash()
  323. return bxor(hash, rshift(hash, 16))
  324. end
  325. function Dictionary:is_empty ()
  326. return self.size == 0
  327. end
  328. function Dictionary:get_bucket_idx (hash)
  329. return band(self.buckets.n - 1, hash) + 1
  330. end
  331. function Dictionary:get_bucket (hash)
  332. return self.buckets[self:get_bucket_idx(hash)]
  333. end
  334. function Dictionary:at (key)
  335. local hash = self:hash(key)
  336. local e = self:get_bucket(hash)
  337. while e do
  338. if e:match(hash, key) then
  339. return e.value
  340. end
  341. e = e.next
  342. end
  343. return nil
  344. end
  345. function Dictionary:contains_key (key)
  346. local hash = self:hash(key)
  347. local e = self:get_bucket(hash)
  348. while e do
  349. if e.match(hash, key) then
  350. return true
  351. end
  352. e = e.next
  353. end
  354. return false
  355. end
  356. function Dictionary:at_put (key, value)
  357. local hash = self:hash(key)
  358. local i = self:get_bucket_idx(hash)
  359. local current = self.buckets[i]
  360. if not current then
  361. self.buckets[i] = self:new_entry(key, value, hash)
  362. self.size = self.size + 1
  363. else
  364. self:insert_bucket_entry(key, value, hash, current)
  365. end
  366. if self.size > self.buckets.n then
  367. self:resize()
  368. end
  369. end
  370. function Dictionary:new_entry (key, value, hash)
  371. return Entry.new(hash, key, value, nil)
  372. end
  373. function Dictionary:insert_bucket_entry (key, value, hash, head)
  374. local current = head
  375. while true do
  376. if current:match(hash, key) then
  377. current.value = value
  378. return
  379. end
  380. if not current.next then
  381. self.size = self.size + 1
  382. current.next = self:new_entry(key, value, hash)
  383. return
  384. end
  385. current = current.next
  386. end
  387. end
  388. function Dictionary:resize ()
  389. local old_storage = self.buckets
  390. self.buckets = alloc_array(old_storage.n * 2)
  391. self:transfer_entries(old_storage)
  392. end
  393. function Dictionary:transfer_entries (old_storage)
  394. local buckets = self.buckets
  395. for i = 1, old_storage.n do
  396. local current = old_storage[i]
  397. if current then
  398. old_storage[i] = nil
  399. if not current.next then
  400. local hash = band(current.hash, buckets.n - 1) + 1
  401. buckets[hash] = current
  402. else
  403. self:split_bucket(old_storage, i, current)
  404. end
  405. end
  406. end
  407. end
  408. function Dictionary:split_bucket (old_storage, i, head)
  409. local lo_head, lo_tail = nil, nil
  410. local hi_head, hi_tail = nil, nil
  411. local current = head
  412. while current do
  413. if band(current.hash, old_storage.n) == 0 then
  414. if not lo_tail then
  415. lo_head = current
  416. else
  417. lo_tail.next = current
  418. end
  419. lo_tail = current
  420. else
  421. if not hi_tail then
  422. hi_head = current
  423. else
  424. hi_tail.next = current
  425. end
  426. hi_tail = current
  427. end
  428. current = current.next
  429. end
  430. if lo_tail then
  431. lo_tail.next = nil
  432. self.buckets[i] = lo_head
  433. end
  434. if hi_tail then
  435. hi_tail.next = nil
  436. self.buckets[i + old_storage.n] = hi_head
  437. end
  438. end
  439. function Dictionary:remove_all ()
  440. self.buckets = alloc_array(self.buckets.n)
  441. self.size = 0
  442. end
  443. function Dictionary:keys ()
  444. local keys = Vector.new(self.size)
  445. local buckets = self.buckets
  446. for i = 1, buckets.n do
  447. local current = buckets[i]
  448. while current do
  449. keys:append(current.key)
  450. current = current.next
  451. end
  452. end
  453. return keys
  454. end
  455. function Dictionary:values ()
  456. local vals = Vector.new(self.size)
  457. local buckets = self.buckets
  458. for i = 1, buckets.n do
  459. local current = buckets[i]
  460. while current do
  461. vals:append(current.value)
  462. current = current.next
  463. end
  464. end
  465. return vals
  466. end
  467. end -- class Dictionary
  468. local IdEntry = {_CLASS = 'IdEntry'} do
  469. setmetatable(IdEntry, {__index = Entry})
  470. function IdEntry.new (hash, key, value, next)
  471. local obj = Entry.new (hash, key, value, next)
  472. return setmetatable(obj, {__index = IdEntry})
  473. end
  474. function IdEntry:match (hash, key)
  475. return self.hash == hash and self.key == key
  476. end
  477. end -- class IdEntry
  478. local IdentityDictionary = {_CLASS = 'IdentityDictionary'} do
  479. setmetatable(IdentityDictionary, {__index = Dictionary})
  480. function IdentityDictionary.new (size)
  481. local obj = Dictionary.new (size)
  482. return setmetatable(obj, {__index = Dictionary})
  483. end
  484. function IdentityDictionary:new_entry (key, value, hash)
  485. return IdEntry.new(hash, key, value, nil)
  486. end
  487. end -- class IdentityDictionary
  488. local Random = {_CLASS = 'Random'} do
  489. function Random.new ()
  490. local obj = {seed = 74755}
  491. return setmetatable(obj, {__index = Random})
  492. end
  493. function Random:next ()
  494. self.seed = band(((self.seed * 1309) + 13849), 65535);
  495. return self.seed;
  496. end
  497. end -- class Random
  498. ---------------------------------
  499. local benchmark = {} do
  500. function benchmark:inner_benchmark_loop (inner_iterations)
  501. for _ = 1, inner_iterations do
  502. if not self:verify_result(self:benchmark()) then
  503. return false
  504. end
  505. end
  506. return true
  507. end
  508. function benchmark:benchmark ()
  509. error 'subclass_responsibility'
  510. end
  511. function benchmark:verify_result ()
  512. error 'subclass_responsibility'
  513. end
  514. end
  515. ---------------------------------
  516. local MIN_X = 0.0
  517. local MIN_Y = 0.0
  518. local MAX_X = 1000.0
  519. local MAX_Y = 1000.0
  520. local MIN_Z = 0.0
  521. local MAX_Z = 10.0
  522. local PROXIMITY_RADIUS = 1.0
  523. local GOOD_VOXEL_SIZE = PROXIMITY_RADIUS * 2.0
  524. local Vector2D = {_CLASS = 'Vector2D'} do
  525. function Vector2D.new (x, y)
  526. local obj = {x = x, y = y}
  527. return setmetatable(obj, {__index = Vector2D})
  528. end
  529. function Vector2D:plus (other)
  530. return Vector2D.new(self.x + other.x, self.y + other.y)
  531. end
  532. function Vector2D:minus (other)
  533. return Vector2D.new(self.x - other.x, self.y - other.y)
  534. end
  535. local function compare_numbers (a, b)
  536. if a == b then
  537. return 0
  538. elseif a < b then
  539. return -1
  540. elseif a > b then
  541. return 1
  542. -- We say that NaN is smaller than non-NaN.
  543. elseif a == a then
  544. return 1
  545. else
  546. return -1
  547. end
  548. end
  549. function Vector2D:compare_to (other)
  550. local result = compare_numbers(self.x, other.x)
  551. if result ~= 0 then
  552. return result
  553. else
  554. return compare_numbers(self.y, other.y)
  555. end
  556. end
  557. end -- class Vector2D
  558. local Vector3D = {_CLASS = 'Vector3D'} do
  559. local sqrt = math.sqrt
  560. function Vector3D.new (x, y, z)
  561. local obj = {x = x, y = y, z = z}
  562. return setmetatable(obj, {__index = Vector3D})
  563. end
  564. function Vector3D:plus (other)
  565. return Vector3D.new(self.x + other.x, self.y + other.y, self.z + other.z)
  566. end
  567. function Vector3D:minus (other)
  568. return Vector3D.new(self.x - other.x, self.y - other.y, self.z - other.z)
  569. end
  570. function Vector3D:dot (other)
  571. return self.x * other.x + self.y * other.y + self.z * other.z
  572. end
  573. function Vector3D:squared_magnitude ()
  574. return self:dot(self)
  575. end
  576. function Vector3D:magnitude ()
  577. return sqrt(self:squared_magnitude())
  578. end
  579. function Vector3D:times (amount)
  580. return Vector3D.new(self.x * amount, self.y * amount, self.z * amount)
  581. end
  582. end -- class Vector3D
  583. local function tree_minimum (x)
  584. local current = x
  585. while current.left do
  586. current = current.left
  587. end
  588. return current
  589. end
  590. local Node = {_CLASS = 'Node'} do
  591. function Node.new (key, value)
  592. local obj = {
  593. key = key,
  594. value = value,
  595. left = nil,
  596. right = nil,
  597. parent = nil,
  598. color = 'red',
  599. }
  600. return setmetatable(obj, {__index = Node})
  601. end
  602. function Node:successor ()
  603. local x = self
  604. if x.right then
  605. return tree_minimum(x.right)
  606. end
  607. local y = x.parent
  608. while y and x == y.right do
  609. x = y
  610. y = y.parent
  611. end
  612. return y
  613. end
  614. end -- class Node
  615. local RbtEntry = {_CLASS = 'RbtEntry'} do
  616. function RbtEntry.new(key, value)
  617. local obj = {
  618. key = key,
  619. value = value
  620. }
  621. return setmetatable(obj, {__index = RbtEntry})
  622. end
  623. end -- class RbtEntry
  624. local InsertResult = {_CLASS = 'InsertResult'} do
  625. function InsertResult.new(is_new_entry, new_node, old_value)
  626. local obj = {
  627. is_new_entry = is_new_entry,
  628. new_node = new_node,
  629. old_value = old_value,
  630. }
  631. return setmetatable(obj, {__index = InsertResult})
  632. end
  633. end -- class InsertResult
  634. local RedBlackTree = {_CLASS = 'RedBlackTree'} do
  635. function RedBlackTree.new ()
  636. local obj = {root = nil}
  637. return setmetatable(obj, {__index = RedBlackTree})
  638. end
  639. function RedBlackTree:put (key, value)
  640. local insertion_result = self:tree_insert(key, value)
  641. if not insertion_result.is_new_entry then
  642. return insertion_result.old_value
  643. end
  644. local x = insertion_result.new_node
  645. while x ~= self.root and x.parent.color == 'red' do
  646. if x.parent == x.parent.parent.left then
  647. local y = x.parent.parent.right
  648. if y and y.color == 'red' then
  649. -- Case 1
  650. x.parent.color = 'black'
  651. y.color = 'black'
  652. x.parent.parent.color = 'red'
  653. x = x.parent.parent
  654. else
  655. if x == x.parent.right then
  656. -- Case 2
  657. x = x.parent
  658. self:left_rotate(x)
  659. end
  660. -- Case 3
  661. x.parent.color = 'black'
  662. x.parent.parent.color = 'red'
  663. self:right_rotate(x.parent.parent)
  664. end
  665. else
  666. -- Same as "then" clause with "right" and "left" exchanged.
  667. local y = x.parent.parent.left
  668. if y and y.color == 'red' then
  669. -- Case 1
  670. x.parent.color = 'black'
  671. y.color = 'black'
  672. x.parent.parent.color = 'red'
  673. x = x.parent.parent
  674. else
  675. if x == x.parent.left then
  676. -- Case 2
  677. x = x.parent
  678. self:right_rotate(x)
  679. end
  680. -- Case 3
  681. x.parent.color = 'black'
  682. x.parent.parent.color = 'red'
  683. self:left_rotate(x.parent.parent)
  684. end
  685. end
  686. end
  687. self.root.color = 'black'
  688. return nil
  689. end
  690. function RedBlackTree:remove (key)
  691. local z = self:find_node(key)
  692. if not z then
  693. return nil
  694. end
  695. -- Y is the node to be unlinked from the tree.
  696. local y
  697. if not z.left or not z.right then
  698. y = z
  699. else
  700. y = z:successor()
  701. end
  702. -- Y is guaranteed to be non-null at this point.
  703. local x
  704. if y.left then
  705. x = y.left
  706. else
  707. x = y.right
  708. end
  709. -- X is the child of y which might potentially replace y
  710. -- in the tree. X might be null at this point.
  711. local x_parent
  712. if x then
  713. x.parent = y.parent
  714. x_parent = x.parent
  715. else
  716. x_parent = y.parent
  717. end
  718. if not y.parent then
  719. self.root = x
  720. else
  721. if y == y.parent.left then
  722. y.parent.left = x
  723. else
  724. y.parent.right = x
  725. end
  726. end
  727. if y ~= z then
  728. if y.color == 'black' then
  729. self:remove_fixup(x, x_parent)
  730. end
  731. y.parent = z.parent
  732. y.color = z.color
  733. y.left = z.left
  734. y.right = z.right
  735. if z.left then
  736. z.left.parent = y
  737. end
  738. if z.right then
  739. z.right.parent = y
  740. end
  741. if z.parent then
  742. if z.parent.left == z then
  743. z.parent.left = y
  744. else
  745. z.parent.right = y
  746. end
  747. else
  748. self.root = y
  749. end
  750. elseif y.color == 'black' then
  751. self:remove_fixup(x, x_parent)
  752. end
  753. return z.value
  754. end
  755. function RedBlackTree:get (key)
  756. local node = self:find_node(key)
  757. if node then
  758. return node.value
  759. end
  760. return nil
  761. end
  762. function RedBlackTree:for_each (fn)
  763. if not self.root then
  764. return
  765. end
  766. local current = tree_minimum(self.root)
  767. while current do
  768. fn(RbtEntry.new(current.key, current.value))
  769. current = current:successor()
  770. end
  771. end
  772. function RedBlackTree:find_node (key)
  773. local current = self.root
  774. while current do
  775. local comparison_result = key:compare_to(current.key)
  776. if comparison_result == 0 then
  777. return current
  778. elseif comparison_result < 0 then
  779. current = current.left
  780. else
  781. current = current.right
  782. end
  783. end
  784. return nil
  785. end
  786. function RedBlackTree:tree_insert (key, value)
  787. local y = nil
  788. local x = self.root
  789. while x do
  790. y = x
  791. local comparison_result = key:compare_to(x.key)
  792. if comparison_result < 0 then
  793. x = x.left
  794. elseif comparison_result > 0 then
  795. x = x.right
  796. else
  797. local old_value = x.value
  798. x.value = value
  799. return InsertResult.new(false, nil, old_value)
  800. end
  801. end
  802. local z = Node.new(key, value)
  803. z.parent = y
  804. if not y then
  805. self.root = z
  806. else
  807. if key:compare_to(y.key) < 0 then
  808. y.left = z
  809. else
  810. y.right = z
  811. end
  812. end
  813. return InsertResult.new(true, z, nil)
  814. end
  815. function RedBlackTree:left_rotate (x)
  816. local y = x.right
  817. -- Turn y's left subtree into x's right subtree.
  818. x.right = y.left
  819. if y.left then
  820. y.left.parent = x
  821. end
  822. -- Link x's parent to y.
  823. y.parent = x.parent
  824. if not x.parent then
  825. self.root = y
  826. else
  827. if x == x.parent.left then
  828. x.parent.left = y
  829. else
  830. x.parent.right = y
  831. end
  832. end
  833. -- Put x on y's left.
  834. y.left = x
  835. x.parent = y
  836. return y
  837. end
  838. function RedBlackTree:right_rotate (y)
  839. local x = y.left
  840. -- Turn x's right subtree into y's left subtree.
  841. y.left = x.right
  842. if x.right then
  843. x.right.parent = y
  844. end
  845. -- Link y's parent to x.
  846. x.parent = y.parent
  847. if not y.parent then
  848. self.root = x
  849. else
  850. if y == y.parent.left then
  851. y.parent.left = x
  852. else
  853. y.parent.right = x
  854. end
  855. end
  856. x.right = y
  857. y.parent = x
  858. return x
  859. end
  860. function RedBlackTree:remove_fixup (x, x_parent)
  861. while x ~= self.root and (not x or x.color == 'black') do
  862. if x == x_parent.left then
  863. -- Note: the text points out that w cannot be null.
  864. -- The reason is not obvious from simply looking at the code;
  865. -- it comes about from the properties of the red-black tree.
  866. local w = x_parent.right
  867. if w.color == 'red' then
  868. -- Case 1
  869. w.color = 'black'
  870. x_parent.color = 'red'
  871. self:left_rotate(x_parent)
  872. w = x_parent.right
  873. end
  874. if (not w.left or w.left.color == 'black') and
  875. (not w.right or w.right.color == 'black') then
  876. -- Case 2
  877. w.color = 'red'
  878. x = x_parent
  879. x_parent = x.parent
  880. else
  881. if not w.right or w.right.color == 'black' then
  882. -- Case 3
  883. w.left.color = 'black'
  884. w.color = 'red'
  885. self:right_rotate(w)
  886. w = x_parent.right
  887. end
  888. -- Case 4
  889. w.color = x_parent.color
  890. x_parent.color = 'black'
  891. if w.right then
  892. w.right.color = 'black'
  893. end
  894. self:left_rotate(x_parent)
  895. x = self.root
  896. x_parent = x.parent
  897. end
  898. else
  899. -- Same as "then" clause with "right" and "left" exchanged.
  900. local w = x_parent.left
  901. if w.color == 'red' then
  902. -- Case 1
  903. w.color = 'black'
  904. x_parent.color = 'red'
  905. self:right_rotate(x_parent)
  906. w = x_parent.left
  907. end
  908. if (not w.right or w.right.color == 'black') and
  909. (not w.left or w.left.color == 'black') then
  910. -- Case 2
  911. w.color = 'red'
  912. x = x_parent
  913. x_parent = x.parent
  914. else
  915. if not w.left or w.left.color == 'black' then
  916. -- Case 3
  917. w.right.color = 'black'
  918. w.color = 'red'
  919. self:left_rotate(w)
  920. w = x_parent.left
  921. end
  922. -- Case 4
  923. w.color = x_parent.color
  924. x_parent.color = 'black'
  925. if w.left then
  926. w.left.color = 'black'
  927. end
  928. self:right_rotate(x_parent)
  929. x = self.root
  930. x_parent = x.parent
  931. end
  932. end
  933. end
  934. if x then
  935. x.color = 'black'
  936. end
  937. end
  938. end -- class RedBlackTree
  939. local CallSign = {_CLASS = 'CallSign'} do
  940. function CallSign.new (value)
  941. local obj = {value = value}
  942. return setmetatable(obj, {__index = CallSign})
  943. end
  944. function CallSign:compare_to (other)
  945. return (self.value == other.value) and 0 or ((self.value < other.value) and -1 or 1)
  946. end
  947. end -- class CallSign
  948. local Collision = {_CLASS = 'Collision'} do
  949. function Collision.new (aircraft_a, aircraft_b, position)
  950. local obj = {
  951. aircraft_a = aircraft_a,
  952. aircraft_b = aircraft_b,
  953. position = position
  954. }
  955. return setmetatable(obj, {__index = Collision})
  956. end
  957. end -- class Collision
  958. local Motion = {_CLASS = 'Motion'} do
  959. local sqrt = math.sqrt
  960. function Motion.new (callsign, pos_one, pos_two)
  961. local obj = {
  962. callsign = callsign,
  963. pos_one = pos_one,
  964. pos_two = pos_two,
  965. }
  966. return setmetatable(obj, {__index = Motion})
  967. end
  968. function Motion:delta ()
  969. return self.pos_two:minus(self.pos_one)
  970. end
  971. function Motion:find_intersection (other)
  972. local init1 = self.pos_one
  973. local init2 = other.pos_one
  974. local vec1 = self:delta()
  975. local vec2 = other:delta()
  976. local radius = PROXIMITY_RADIUS
  977. -- this test is not geometrical 3-d intersection test,
  978. -- it takes the fact that the aircraft move
  979. -- into account; so it is more like a 4d test
  980. -- (it assumes that both of the aircraft have a constant speed
  981. -- over the tested interval)
  982. -- we thus have two points,
  983. -- each of them moving on its line segment at constant speed;
  984. -- we are looking for times when the distance between
  985. -- these two points is smaller than r
  986. -- vec1 is vector of aircraft 1
  987. -- vec2 is vector of aircraft 2
  988. -- a = (V2 - V1)^T * (V2 - V1)
  989. local a = vec2:minus(vec1):squared_magnitude()
  990. if a ~= 0.0 then
  991. -- we are first looking for instances
  992. -- of time when the planes are exactly r from each other
  993. -- at least one plane is moving;
  994. -- if the planes are moving in parallel, they do not have constant speed
  995. -- if the planes are moving in parallel, then
  996. -- if the faster starts behind the slower,
  997. -- we can have 2, 1, or 0 solutions
  998. -- if the faster plane starts in front of the slower,
  999. -- we can have 0 or 1 solutions
  1000. -- if the planes are not moving in parallel, then
  1001. -- point P1 = I1 + vV1
  1002. -- point P2 = I2 + vV2
  1003. -- - looking for v, such that dist(P1,P2) = || P1 - P2 || = r
  1004. -- it follows that || P1 - P2 || = sqrt( < P1-P2, P1-P2 > )
  1005. -- 0 = -r^2 + < P1 - P2, P1 - P2 >
  1006. -- from properties of dot product
  1007. -- 0 = -r^2 + <I1-I2,I1-I2> + v * 2<I1-I2, V1-V2> + v^2 *<V1-V2,V1-V2>
  1008. -- so we calculate a, b, c - and solve the quadratic equation
  1009. -- 0 = c + bv + av^2
  1010. -- b = 2 * <I1-I2, V1-V2>
  1011. local b = 2.0 * init1:minus(init2):dot(vec1:minus(vec2))
  1012. -- c = -r^2 + (I2 - I1)^T * (I2 - I1)
  1013. local c = -radius * radius + init2:minus(init1):squared_magnitude()
  1014. local discr = b * b - 4.0 * a * c
  1015. if discr < 0.0 then
  1016. return nil
  1017. end
  1018. local v1 = (-b - sqrt(discr)) / (2.0 * a)
  1019. local v2 = (-b + sqrt(discr)) / (2.0 * a)
  1020. if v1 <= v2 and ((v1 <= 1.0 and 1.0 <= v2) or
  1021. (v1 <= 0.0 and 0.0 <= v2) or
  1022. (0.0 <= v1 and v2 <= 1.0)) then
  1023. -- Pick a good "time" at which to report the collision.
  1024. local v
  1025. if v1 <= 0.0 then
  1026. -- The collision started before this frame.
  1027. -- Report it at the start of the frame.
  1028. v = 0.0
  1029. else
  1030. -- The collision started during this frame. Report it at that moment.
  1031. v = v1
  1032. end
  1033. local result1 = init1:plus(vec1:times(v))
  1034. local result2 = init2:plus(vec2:times(v))
  1035. local result = result1:plus(result2):times(0.5)
  1036. if result.x >= MIN_X and
  1037. result.x <= MAX_X and
  1038. result.y >= MIN_Y and
  1039. result.y <= MAX_Y and
  1040. result.z >= MIN_Z and
  1041. result.z <= MAX_Z then
  1042. return result
  1043. end
  1044. end
  1045. return nil
  1046. end
  1047. -- the planes have the same speeds and are moving in parallel
  1048. -- (or they are not moving at all)
  1049. -- they thus have the same distance all the time;
  1050. -- we calculate it from the initial point
  1051. -- dist = || i2 - i1 || = sqrt( ( i2 - i1 )^T * ( i2 - i1 ) )
  1052. local dist = init2:minus(init1):magnitude()
  1053. if dist <= radius then
  1054. return init1:plus(init2):times(0.5)
  1055. end
  1056. return nil
  1057. end
  1058. end -- class Motion
  1059. local CollisionDetector = {_CLASS = 'CollisionDetector'} do
  1060. local floor = math.floor
  1061. local HORIZONTAL = Vector2D.new(GOOD_VOXEL_SIZE, 0.0)
  1062. local VERTICAL = Vector2D.new(0.0, GOOD_VOXEL_SIZE)
  1063. function CollisionDetector.new ()
  1064. local obj = {state = RedBlackTree.new()}
  1065. return setmetatable(obj, {__index = CollisionDetector})
  1066. end
  1067. function CollisionDetector:handle_new_frame (frame)
  1068. local motions = Vector.new()
  1069. local seen = RedBlackTree.new()
  1070. frame:each(function (aircraft)
  1071. local old_position = self.state:put(aircraft.callsign, aircraft.position)
  1072. local new_position = aircraft.position
  1073. seen:put(aircraft.callsign, true)
  1074. if not old_position then
  1075. -- Treat newly introduced aircraft as if they were stationary.
  1076. old_position = new_position
  1077. end
  1078. motions:append(Motion.new(aircraft.callsign, old_position, new_position))
  1079. end)
  1080. -- Remove aircraft that are no longer present.
  1081. local to_remove = Vector.new()
  1082. self.state:for_each(function (e)
  1083. if not seen:get(e.key) then
  1084. to_remove:append(e.key)
  1085. end
  1086. end)
  1087. to_remove:each(function (e)
  1088. self.state:remove(e)
  1089. end)
  1090. local all_reduced = self:reduce_collision_set(motions)
  1091. local collisions = Vector.new()
  1092. all_reduced:each(function (reduced)
  1093. for i = 1, reduced:size() do
  1094. local motion1 = reduced:at(i)
  1095. for j = i + 1, reduced:size() do
  1096. local motion2 = reduced:at(j)
  1097. local collision = motion1:find_intersection(motion2)
  1098. if collision then
  1099. collisions:append(Collision.new(motion1.callsign,
  1100. motion2.callsign,
  1101. collision))
  1102. end
  1103. end
  1104. end
  1105. end)
  1106. return collisions
  1107. end
  1108. function CollisionDetector:is_in_voxel (voxel, motion)
  1109. if voxel.x > MAX_X or
  1110. voxel.x < MIN_X or
  1111. voxel.y > MAX_Y or
  1112. voxel.y < MIN_Y then
  1113. return false
  1114. end
  1115. local init = motion.pos_one
  1116. local fin = motion.pos_two
  1117. local v_s = GOOD_VOXEL_SIZE
  1118. local r = PROXIMITY_RADIUS / 2.0
  1119. local v_x = voxel.x
  1120. local x0 = init.x
  1121. local xv = fin.x - init.x
  1122. local v_y = voxel.y
  1123. local y0 = init.y
  1124. local yv = fin.y - init.y
  1125. local low_x = (v_x - r - x0) / xv
  1126. local high_x = (v_x + v_s + r - x0) / xv
  1127. if xv < 0.0 then
  1128. local tmp = low_x
  1129. low_x = high_x
  1130. high_x = tmp
  1131. end
  1132. local low_y = (v_y - r - y0) / yv
  1133. local high_y = (v_y + v_s + r - y0) / yv
  1134. if yv < 0.0 then
  1135. local tmp = low_y
  1136. low_y = high_y
  1137. high_y = tmp
  1138. end
  1139. return (((xv == 0.0 and v_x <= x0 + r and x0 - r <= v_x + v_s) or -- no motion in x
  1140. (low_x <= 1.0 and 1.0 <= high_x) or (low_x <= 0.0 and 0.0 <= high_x) or
  1141. (0.0 <= low_x and high_x <= 1.0)) and
  1142. ((yv == 0.0 and v_y <= y0 + r and y0 - r <= v_y + v_s) or -- no motion in y
  1143. ((low_y <= 1.0 and 1.0 <= high_y) or (low_y <= 0.0 and 0.0 <= high_y) or
  1144. (0.0 <= low_y and high_y <= 1.0))) and
  1145. (xv == 0.0 or yv == 0.0 or -- no motion in x or y or both
  1146. (low_y <= high_x and high_x <= high_y) or
  1147. (low_y <= low_x and low_x <= high_y) or
  1148. (low_x <= low_y and high_y <= high_x)))
  1149. end
  1150. function CollisionDetector:put_into_map (voxel_map, voxel, motion)
  1151. local array = voxel_map:get(voxel)
  1152. if not array then
  1153. array = Vector.new()
  1154. voxel_map:put(voxel, array)
  1155. end
  1156. array:append(motion)
  1157. end
  1158. function CollisionDetector:recurse (voxel_map, seen, next_voxel, motion)
  1159. if not self:is_in_voxel(next_voxel, motion) then
  1160. return
  1161. end
  1162. if seen:put(next_voxel, true) then
  1163. return
  1164. end
  1165. self:put_into_map(voxel_map, next_voxel, motion)
  1166. self:recurse(voxel_map, seen, next_voxel:minus(HORIZONTAL), motion)
  1167. self:recurse(voxel_map, seen, next_voxel:plus(HORIZONTAL), motion)
  1168. self:recurse(voxel_map, seen, next_voxel:minus(VERTICAL), motion)
  1169. self:recurse(voxel_map, seen, next_voxel:plus(VERTICAL), motion)
  1170. self:recurse(voxel_map, seen, next_voxel:minus(HORIZONTAL):minus(VERTICAL), motion)
  1171. self:recurse(voxel_map, seen, next_voxel:minus(HORIZONTAL):plus(VERTICAL), motion)
  1172. self:recurse(voxel_map, seen, next_voxel:plus(HORIZONTAL):minus(VERTICAL), motion)
  1173. self:recurse(voxel_map, seen, next_voxel:plus(HORIZONTAL):plus(VERTICAL), motion)
  1174. end
  1175. function CollisionDetector:reduce_collision_set (motions)
  1176. local voxel_map = RedBlackTree.new()
  1177. motions:each(function (motion)
  1178. self:draw_motion_on_voxel_map(voxel_map, motion)
  1179. end)
  1180. local result = Vector.new()
  1181. voxel_map:for_each(function (e)
  1182. if e.value:size() > 1 then
  1183. result:append(e.value)
  1184. end
  1185. end)
  1186. return result
  1187. end
  1188. function CollisionDetector:voxel_hash (position)
  1189. local x_div = floor(position.x / GOOD_VOXEL_SIZE)
  1190. local y_div = floor(position.y / GOOD_VOXEL_SIZE)
  1191. local x = GOOD_VOXEL_SIZE * x_div
  1192. local y = GOOD_VOXEL_SIZE * y_div
  1193. if position.x < 0 then
  1194. x = x - GOOD_VOXEL_SIZE
  1195. end
  1196. if position.y < 0 then
  1197. y = y - GOOD_VOXEL_SIZE
  1198. end
  1199. return Vector2D.new(x, y)
  1200. end
  1201. function CollisionDetector:draw_motion_on_voxel_map (voxel_map, motion)
  1202. local seen = RedBlackTree.new()
  1203. return self:recurse(voxel_map, seen, self:voxel_hash(motion.pos_one), motion)
  1204. end
  1205. end -- class CollisionDetector
  1206. local Aircraft = {_CLASS = 'Aircraft'} do
  1207. function Aircraft.new (callsign, position)
  1208. local obj = {
  1209. callsign = callsign,
  1210. position = position
  1211. }
  1212. return setmetatable(obj, {__index = Aircraft})
  1213. end
  1214. end -- class Collision
  1215. local Simulator = {_CLASS = 'Simulator'} do
  1216. local cos = math.cos
  1217. local sin = math.sin
  1218. function Simulator.new (num_aircrafts)
  1219. local aircraft = Vector.new()
  1220. for i = 1, num_aircrafts do
  1221. aircraft:append(CallSign.new(i))
  1222. end
  1223. local obj = {aircraft = aircraft}
  1224. return setmetatable(obj, {__index = Simulator})
  1225. end
  1226. function Simulator:simulate (time)
  1227. local frame = Vector.new()
  1228. for i = 1, self.aircraft:size() - 1, 2 do
  1229. frame:append(Aircraft.new(self.aircraft:at(i),
  1230. Vector3D.new(time,
  1231. cos(time) * 2 + (i - 1) * 3,
  1232. 10.0)))
  1233. frame:append(Aircraft.new(self.aircraft:at(i + 1),
  1234. Vector3D.new(time,
  1235. sin(time) * 2 + (i - 1) * 3,
  1236. 10.0)))
  1237. end
  1238. return frame
  1239. end
  1240. end -- class Simulator
  1241. local cd = {} do
  1242. setmetatable(cd, {__index = benchmark})
  1243. function cd:benchmark (num_aircrafts)
  1244. local num_frames = 200
  1245. local simulator = Simulator.new(num_aircrafts)
  1246. local detector = CollisionDetector.new()
  1247. local actual_collisions = 0
  1248. for i = 0, num_frames - 1 do
  1249. local time = i / 10.0
  1250. local collisions = detector:handle_new_frame(simulator:simulate(time))
  1251. actual_collisions = actual_collisions + collisions:size()
  1252. end
  1253. return actual_collisions
  1254. end
  1255. function cd:inner_benchmark_loop (inner_iterations)
  1256. return self:verify_result(self:benchmark(inner_iterations), inner_iterations)
  1257. end
  1258. function cd:verify_result (actual_collisions, num_aircrafts)
  1259. if num_aircrafts == 1000 then
  1260. return actual_collisions == 14484
  1261. elseif num_aircrafts == 500 then
  1262. return actual_collisions == 14484
  1263. elseif num_aircrafts == 250 then
  1264. return actual_collisions == 10830
  1265. elseif num_aircrafts == 200 then
  1266. return actual_collisions == 8655
  1267. elseif num_aircrafts == 100 then
  1268. return actual_collisions == 4305
  1269. elseif num_aircrafts == 10 then
  1270. return actual_collisions == 390
  1271. elseif num_aircrafts == 2 then
  1272. return actual_collisions == 42
  1273. else
  1274. print(('No verification result for %d found'):format(num_aircrafts))
  1275. print(('Result is: %d'):format(actual_collisions))
  1276. return false
  1277. end
  1278. end
  1279. end -- object cd
  1280. return function(N)
  1281. cd:inner_benchmark_loop(N)
  1282. end