Od pohledu dobrý, aneb jak najít skoro stejné obrázky mezi dvěma miliony souborů za méně než deset minut

A tenhle znáte: „Proč je Pentium rychlejší než 486? 486 počítá, ale Pentium jen odhaduje.“ Přestože jde o nejapný vtip o chybě v prvních Pentiích z počítačového pravěku, má v sobě zrnko pravdy: Odhad může být mnohem rychlejší než přesný výpočet.

Odhad, pokud chci, aby byl k něčemu dobrý, není vůbec jednoduchou záležitostí. Nestačí jen vyprodukovat nějaké hausnumero, ale přibližnou odpověď, která má jasně definovanou chybu. Mezi algoritmy, které dělají právě tohle patří kromě klasických velkých jmen jako Bloomfilter, HyperLogLog, Count-min Skech, MinHash a SimHash, patří také perceptuální hashe určené k detekci téměř shodných obrázků.

Všechny zmíněné algoritmy se spoléhají na hashovací funkce, ale některé (MinHash, SimHash a perceptuální hashe o kterých to dneska celé je) jinak než je běžně známé.

Klasické kryptografické hashe, jako je MD5, SHA1 a jejich kamarádi, jsou navrženy tak, že změna jednoho bitu vstupu v ideálním případě znamená změnu poloviny bitů výstupu. Nekryptografické hashe, běžně používané v programování, nemají tuto vlastnost. Typicky jde aspoň o univerzální hash nebo o nějaký k-nezávislý hash. Otisk je přesto výsledkem divokých permutací vstupních dat a neměl by zachovávat jejich korelace. Perceptuální hash je přesně opačný – jde o takový otisk obrázku, který zachovává korelace vstupních dat – když jsou si obrázky na vstupu podobné, budou si s největší pravděpodobností podobné i jejich otisky a míra podobnosti obrázků odpovídá podobnosti hashů.

Algoritmus je jednoduchý:

  1. Zmenším obrázek například na 32×32 pixelů (dělám to kvůli zrychlení následujících výpočtů a ne pro redukci vysokých frekvencí)
  2. Převedu na odstíny šedi (opět kvůli zrychlení následujících kro­ků)
  3. Spočítám 32×32 diskrétní kosinovou transformaci (DCT)
  4. Zredukuji DCT a nechám si jen prvních 8×8 hodnot (kromě té úplně první, která by rozhodila průměr). Tyto hodnoty reprezentují nejnižší frekvence obrázku, tedy ty, které popisují hrubou strukturu.
  5. Spočítám průměrnou hodnotu těchto frekvencí
  6. Vytvořím finální hash o délce 64 bitů – daný bit má hodnotu 1, když je korespondující frekvence vyšší než průměr, 0, pokud je menší

Takto získám perceptuální hash obrázku, kterému nevadí změny barevnosti, světlosti (kvůli tomu, že zaznamenávám rozdíl proti průměru) a ani určité úpravy obrázku samotného. Nejnižší frekvence zaznamenávají nejhrubší kontury a rozdělení světel a stínů a jsou proto odolné vůči změnám v detailech.

Když pak chci porovnat, jak jsou si dva obrázky podobné, spočítám Hammingovu vzdálenost mezi jejich perceptuálními hashi. Výsledkem je počet bitů, ve kterých se oba hashe liší. Když se neliší v žádném bitu, jde o dva vizuálně téměř identické obrázky. Možná byl jednou uložen jako jpg a podruhé jako png. Když se liší jen pár bitů (±5) jde nejspíš o jeden a tentýž obrázek s drobnými úpravami, když se liší ve více bitech, jde s největší pravděpodobností o rozdílné obrázky.

Zkusmo jsem si perceptuální hashování napsal v jazyce Go (respektive bezostyšně ukradl přepsal tohle z Javy). Kompletní kód je k vidění v samostatném gistu.

func ImagePHash(img image.Image) uint64 {
  c := make([]float64, 32)
  for i := 1; i < 32; i++ {
    c[i] = 1
  }
  c[0] = 1 / math.Sqrt(2.0)

  // Reduce size and reduce color
  img = grayscale(resize.Resize(32, 32, img, resize.Bilinear))

  vals := make([][]float64, 32)
  for i := range vals {
    vals[i] = make([]float64, 32)
  }

  bounds := img.Bounds()
  for x := 0; x < bounds.Max.X; x++ {
    for y := 0; y < bounds.Max.Y; y++ {
      _, _, b, _ := img.At(x,y).RGBA()
      vals[x][y] = float64(b)
    }
  }

  /* Compute the DCT */
  dctVals := applyDCT(vals, c)

  // Reduce the DCT and compute the average value (excluding the first term)
  total := 0.0

  for x := 0; x < 8; x++ {
    for y := 0; y < 8; y++ {
      total += dctVals[x][y]
    }
  }
  total -= dctVals[0][0]

  avg := total / float64((8 * 8) - 1)

  // Further reduce the DCT and create 64 bit hash
  hash := uint64(0)

  for x := 0; x < 8; x++ {
    for y := 0; y < 8; y++ {
      if !(x == 0 && y == 0) {
        if dctVals[x][y] > avg {
          hash = hash | (1 << (63-uint64(x*8+y)))
        }
      }
    }
  }

  return hash
}

func applyDCT(f [][]float64, c []float64) [][]float64 {
  F := make([][]float64, 32)
  for i := range F {
    F[i] = make([]float64, 32)
  }

  cosines := [32][32]float64 {}
  for a := 0; a < 32 ; a++ {
    for b := 0; b < 32 ; b++ {
      cosines[a][b] = math.Cos((float64(2*a+1) / float64(2*32)) * float64(b) * math.Pi)
    }
  }

  for u := 0; u < 32; u++ {
    for v := 0; v < 32; v++ {
      sum := 0.0
      for i := 0; i < 32; i++ {
        for j := 0; j < 32; j++ {
          sum += cosines[i][u] * cosines[j][v] * f[i][j]
        }
      }
      sum *= (c[u]*c[v] / 4.0)
      F[u][v] = sum
    }
  }
  return F
}

Jak je vidět, jde o pár řádků kódu a výsledky jsou překvapivě přesné.

(nebo, nebo, nebo)

Přesto nejde o neomylný algoritmus. Kód si pohrává s pravděpodobnostmi a někdy, když se plete, bývají výsledky veselé. Zvlášť, když zvýším toleranci podobnosti.

Teď už zbývá jen vyřešit, jak to dělat rychle, když chci projít v titulku zmíněné dva miliony obrázků pod deset minut.

Samotné porovnání p-hashů se může zdát jako dostatečně rychlé. V ideálním případě jde o pouhé 2 instrukce: xor a popcnt + něco málo okolo.1 Přesto, pokud chci porovnat všechny hashe se všemi ostatními, musím udělat kvadratické množství práce a velké číslo se umocněním na druhou bohužel nezmenší.

Naštěstí můžu všechno znatelně urychlit pomocí postupu známého jako locality-sensitive hashing (LSH) – další z metod, která si hraje s pravděpodobnostmi a obětuje něco málo přesnosti ve prospěch masivního zrychlení.

LSH využívá hashe zachovávající lokalitu (podobné věci mají podobné hashe) a snaží se položky rozdělit do skupin, kdy v každé skupině jsou položky, které si můžou být podobné, ale nejsou tam ty, které jsou očividně rozdílné.

V tomto případě zkonstruuji LSH pro Hammingovu vzdálenost tak, že 64 bitů hashe rozsekám na 8 pruhů (band) po 8 bitech a každý hash přiřadím do skupiny (bucket) určenou bity každého pruhu.

To znamená, že když prvních osm bitů hashe má hodnotu 0xFF, v prvním pruhu ho přiřadím do bucketu 0xFF, když má druhých 8 bitů 0x0A, v druhém pruhu spadne do skupiny 0x0A a tak dále.

Celé to funguje na jednoduchém předpokladu: Pokud jsou si dva hashe velice podobné, s největší pravděpodobností se shodují v bitech jednoho nebo více pruhů.2 Když chci najít všechny hashe podobné hashi A, vezmu jeho prvních 8 bitů, podívám se do odpovídajícího bucketu prvního pruhu a spočítám podobnost jen s těmito kandidáty, pak vezmu dalších osm bitů a spočítám podobnost s těmito kandidáty. To udělám jednou pro každý pruh a mám výsledné podobné hashe. Protože jde o pravděpodobnostní metodu, může se stát, že některé skutečně podobné hashe to přeskočí, protože se vždy o jeden bit lišily ve všech pruzích a proto skončily v odlišných bucketech a vůbec nebyly vybráni jako kandidáti. Pokud se zajímám o velice podobné hashe, je tato šance velice malá. Pokud chci v takto vytvořeném LSH najít hashe lišící se v méně než 8 bitech, je garantováno, že všechny dostanu, protože se aspoň v jednom pruhu musí shodovat.

S LSH, kdy mám 8 pruhů po 8 bitech, mi musím udělat jen 8/256 práce – 32× méně, protože se musím podívat do 8 bucketů (jeden v každém pruhu) a spočítat podobnosti ve všemi kandidáty, kterých je 1/256. Pokud to nestačí, můžu tyto dva parametry zvolit podle toho, jak je třeba: když chci víc přesnosti zvolím menší šířku pruhu (12×5 bitů), pokud chci více rychlosti, zvolím širší pruhy (6×10 bitů).

No a s tímto „zlepšovákem“ dokážu bez větších problémů porovnat ony dva miliony obrázků pod 10 minut.

pattern := "/path/to/images/*"
fs, _ := filepath.Glob(pattern)

files  := make([]string, 0, len(fs))
hashes := make([]uint64, 0, len(fs))

for _, fullPath := range fs {

  reader, err := os.Open(fullPath)
  if err != nil { continue }

  img, _, err := image.Decode(reader)
  if err != nil { continue }

  hash := ImagePHash(img)
  files  = append(files, fullPath)
  hashes = append(hashes, hash)

  reader.Close()
}


// create LSH map

const Bands = 8
const BandBits = 8
const BandMask = (1 << BandBits) - 1

bands := [Bands][BandMask+1][]IdxHashPair {}

for idx, hash := range hashes {
  for b := 0 ; b < Bands ; b++ {
    h := hash << uint32(b*BandBits) & BandMask
    bands[b][h] = append(bands[b][h], IdxHashPair { idx, hash })
  }
}


// compute similarities

for idx, hash := range hashes {
  for b := 0 ; b < Bands ; b++ {
    h := hash << uint32(b*BandBits) & BandMask

    for _, candidate := range bands[b][h] {
      if idx <= candidate.Idx { continue }

      diff := differentBits(hash, candidate.Hash)
      if diff < 6  {
        similarity := 1.0 - (float64(diff) / float64(64))
        fmt.Printf("%s %s %d\n", files[idx], files[candidate.Idx], similarity)
      }
    }
  }
}

To je všechno. Záhada vyřešena.


Protože tuto techniku používám na několika místech (například pro detekci podobných článků na devblozích, tady, na téhle maličkosti a na mnoha dalších místech), tak jsem si proto napsal vlastní knihovnu nazvanou prozaicky sketches, která obsahuje implementaci PHashe, skečů, LSH, ALS, FunkSVD nebo DBSCAN clusterování a dalších nesouvisejících vě­cí.


Dále k tématu:


Pozn:

  1. V Go se bohužel k popcnt bez tanečků v assembleru nedá dostat a tak je třeba si napsat utilitu, která dělá to samé, ale potřebuje o něco víc instrukcí.
  2. Pravděpodobnost, že LSH odhalí dva kandidáty je 1 – (1 – pn)b, kde p je podobnost mezi dvěma kandidáty, n šířka pruhu a b je počet pruhů. Pomocí tohoto vztahu můžu najít hodnoty n a p, které budou fungovat pro požadovaný práh podobnosti.

Flattr this!

Posted in Algo, DS | Leave a comment

Sketches – slajdy

Sketches – slajdy

Flattr this!

Posted in Uncategorized | Leave a comment

YOLO tree

Stará internetová moudrost říká „you only live once“. Život je příliš krátký na jistoty a někdy to prostě musíme risknout. Jako třeba během konstrukce binárního vyhledávacího stromu (BST).

BST je super věc, všichni ho známe, milujeme a napíšeme si jeden nebo dva ještě před snídaní. Ale i přesto má určité vlastnosti, které, i když se zdají naprosto samozřejmé, nemusí zrovna být to, co by člověk chtěl. Konstrukce BST je striktně sekvenční a nedá se (přímo1) naivně paralelizovat. To znamená, že pokud chci nejmenší element nějaké sekvence, musím nejdřív postavit celý BST v čase O(n log(n)) a teprve potom z něj vyhmátnout minimum. V tom okamžiku je ale už pozdě, protože jsem udělal všechnu práci nutnou k seřazení sekvence. V ideálním případě bych udělal jen nezbytně nutné množství práce a konstrukci zbytku bych líně odložil na později.

Naproti tomu je quicksort (QS) nebo radix sort (RS) představuje potenciálně líné algoritmy. Když chci n-tý nejmenší prvek, tak si ho v jednoduše quickselectnu a cestou udělám lineární množství práce vstříc plně seřazenému stupu. Stejně tak jsou všechny QS/RS naivně paralelizovatelé. Jako všechny algoritmy, které vstup rekurzivně rozdělí a pak jednotlivé části nezávisle zpracovávají, se u nich dá triviálně dosáhnout logaritmického paralelního zrychlení.

BST se od QS/RS odlišuje tím, že výsledkem není jen seřazená sekvence, ale dynamická datová struktura, do které můžu přidávat nebo z ní odebírat prvky.

Když chci to nejlepší z obou světů, musím na to jít jako Abe a Aaron z Primeru:

They took from their surroundings what was needed and made of it something more.

Reifikací quicksortu můžu dostat BST, který je dynamický a vyvážený a jeho konstrukce je líná a naivně paralelizovatelná.

Jako při každé reifikaci začnu s nehmotným procesem a zhmotním ho do datové struktury. Když v QS vybírám pivot, reifikuji ho jako uzel stromu s tím, že všechny menší prvky přijdou doleva a ty větší doprava. Takhle postupuji rekurzivně k listům a, když si to moje srdce žádá, strom průběžně vyvažuji.

Strom je tvořen dvěma druhy uzlů – nerealizovanými a realizovanými. Ty nerealizované jsou jednoduché seznamy prvků, realizované uzly jsou pak tvořeny pivotem a dvěma podstromy.

Začnu s nerealizovaným seznamem:

[6, 8, 3, 1, 4, 2, 9, 5, 0, 7]

Vyberu pivot, udělám z něj jeden realizovaný uzel, rozdělím seznam na dva nerealizované seznamy a ty strčím jako podstromy. V tomto případě výběr probíhá v Hoarově stylu první dobrá, ale je možné použít i komplikovanější metody jako best-of-3 nebo quasi best-of-9.

                   (6)
                  /   \
[3, 1, 4, 2, 5, 0]     [8, 9, 7]

(Nerealizované seznamy uvádím vždy v hranatých závorkách. Do kulatých píšu pivot realizovaného uz­lu.)

Takhle můžu pokračovat, dokud budu chtít. Můžu třeba realizovat jen cestu k nejmenšímu prvku v průměrném lineárním čase:

            (6)
           /   \
        (3)     [8, 9, 7]
       /   \
    (1)     [4, 5]
   /   \
(0)     [2]

Jak je vidět, velká část stromu vůbec nebyla vyhodnocena.

Stejně tak můžu strom postavit najednou (podstromy můžu potenciálně stavět paralelně) rekurzivní aplikací tohoto pravidla.

Tato QS inspirovaná konstrukce je úzce spjatá s BST. Knuth o tom psal v The Art of Computer Programming: Když jako pivot používám vždycky první element, je výsledný strom stejný jako když stavím binární strom bez vyvažování. Ve výsledku provedu stejná porovnání, jen v jiném pořadí.

Pokud strom nesplňuje moje požadavky na vyváženost, můžu ho průběžně balancovat.

V mém případě může být levý podstrom příliš velký:

                   (6)
                  /   \
[3, 1, 4, 2, 5, 0]     [8, 9, 7]

Rotace si v tomto případě vynutí částečnou realizaci podstromu.

              (6)
             /   \
          (3)     [8, 9, 7]
         /   \
[1, 0, 2]     (4)
                 \
                  [5]

Po rotaci to bude vypadat takhle:

              (4)
             /   \
          (3)     (6)
         /       /   \
[1, 0, 2]     [5]     [8, 9, 7]

Musel jsme realizovat jednu úroveň stromu. Tento proces může pokračovat rekurzivně až k listům v průměrném lineárním čase (a quicksortovské O(n2) v nejhorším případě).

K vyvažování se nehodí techniky z red-black nebo AVL stromů. Ty potřebují označit každý uzel (buď barvou nebo hloubkou) a to nutně vede k realizaci celého podstromu.

Když jsem o tom přemýšlel, nějakou dobu jsem nevěděl, jak tenhle problém vyřešit. Pak jsem ale v paperu Implementing sets efficiently in a functional language narazil na size balanced trees – tedy stromy vyvažované podle velikosti. Velikost podstromu je lokální vlastnost, kterou znám nehledě na to, zdali je podstrom realizovaný nebo ne. Balancování začne jen, když je jeden podstrom k-krát větší než ten druhý.

Vyvažování ale v principu nemusí být nutné, protože když vybírám pivot náhodně, v průměrném případě je výsledný strom vybalancovaný s průměrnou hloubkou 1.386 log(n). Takže místo balancování bych mohl jen randomizovat pořadí vstupu.

Quicksort styl konstrukce vede k lenosti. Teď jak na dynamické vlastnosti? Chci mít možnost do stromu přidávat a odebírat prvky a zachovat při tom jeho lenost a asymptoty.

Vkládání je triviální. Místo toho, aby byl uzel buď realizovaný nebo ne, může být realizovaný se seznam nerealizovaných změn.

Když například do tohohle stromu

            (6)
           /   \
        (3)     [8, 9, 7]
       /   \
    (1)     [4, 5]
   /   \
(0)     [2]

vložím [ 1.47 4.7 47 ], dostanu

            (6)-[1.47, 4.7, 47]
           /   \
        (3)     [8, 9, 7]
       /   \
    (1)     [4, 5]
   /   \
(0)     [2]

Když budu opět chtít získat minimum, všechny nerealizované seznamy, na které narazím cestou, zrealizuji:

            (6)
           /   \
        (3)     [8, 9, 7, 47]
       /   \
    (1)     [4, 5, 4.7]
   /   \
(0)     [2, 1.47]

V tomto případě jsem jen přidával na konec nerealizovaných seznamů.

Mazání můžu vyřešit dvěma způsoby. Buď striktně nebo líně. Striktní varianta najde mazaný prvek ve stromu a hned ho odstraní. Cestou nutně realizuje všechny uzly (to by mělo být v nejhorším případě O(n) práce). Líná varianta do nerealizovaného seznamu jen přidá tombstone. Když se tento objekt eventuálně dostane do listu, odstaní daný element. Tohle řešení je líné, ale zase může vést k hromadění tombstone objektů v nikdy nerealizovaných částech stromu. To se dá vyřešit tak, že každý podstrom nejen sleduje svou velikost, ale také počet náhrobků.

A to by bylo všechno.


Do tohoto bodu jsem se dostal jenom tak, že jsem přemýšlel, kam až se dá výchozí myšlenka dotáhnout. Pak jsem si ale najednou uvědomil, že jsem už něco podobného viděl.

Fraktální stromy, používané v TokuDB, mají velice podobnou strukturu. Druhá verze fraktálních stromů jsou B-stromy, kde je ke každému vnitřnímu uzlu připnutý buffer akumulující změny. Ty jsou eventuálně (např. když se buffer naplní), propagovány o úroveň níž do bufferů v příslušných uzlech nebo přímo do listů s datovými položkami.

Fraktální stromy jsou podobné mým YOLO stromům, i když k této podobě došly z jiného úhlu. TokuDB se snažilo omezit náhodné IO operace a nahradit je sekvenčním zápisem. Když se buffer naplní, jeho čtení a zápis do nižších úrovní je striktně sekvenční. Já jsem se v případě YOLO stromů zajímal jen o líné zpracování a sledoval ho až do extrému. V obou případech to došlo k podobným závěrům.

Pokud jsou tedy YOLO/fraktal tree reifikacemi quicksortu (resp. multi-way quicksortu), musí existovat podobná konstrukce, která reifikuje merge sort. Ta skutečně existuje a co víc – používá se ve světě databází.

Jde o LSMT – Log Structured Merge Trees v kombinaci s tzv. fractional cascasing – reifikovaný merge sort, který je líný a má ty správné asymptoty jako například O(log(n)) hledání.

Shodou náhod na LSMT byla založena první verze fraktálních stromů z TokuDB.


Dále k tématu:


  1. Můžu paralelně postavit několik BST a pak je naivně sloučit s maximálně logaritmickým zrychlením.

Flattr this!

Posted in Algo | Leave a comment

diff a stromy

Minule jsem tu psal, jak dělat diffy a jak se tato znalost dá využít pro kompresi dat. Teď bych chtěl jít o krok dál k diffování libovolných stromů.

Už tehdy jsem se zmínil, že jsem o tom chvíli přemýšlel, ale jako vždycky jsem se nedostal příliš daleko. V hlavě jsem měl jen řešení několika triviálních případů, ale žádné principy. Pak ale moje obsesivní povaha převzala kontrolu a tendence všechno dotáhnout tak daleko, jak jen to bude možné. Jako vždycky jsem potřeboval algoritmus.

Když se dostanu do problémů, měl bych vždycky se vždycky zeptat WWEKD – What would Edward Kmett do? On má řešení mnoha problémů, tak proč ne zrovna toho současného?

Když jsem se rozhlížel po netu, zjistil jsem, že existují nějaké knihovny, jako třeba diffson, které počítají diffy JSON stromů. Ty ale neumí všechno. diffson dělá diffy jen úroveň po úrovni. Zkontroluje kořen jednoho stromu s druhým a když se liší, začne spolu porovnávat vnořené uzly. Neumí uspokojivě spočítat diff pro rozsáhlejší změny struktury stromu.

Například mám strom

        (5)
       /   \
    (2)     [6]
   /   \
[1]     [3]

a chci zjistit jeho diff se stromem, který jen jinak zcela schodný, ale má jen jednu úroveň navíc.

    (4)
   /   \
[0]     (5)
       /   \
    (2)     [6]
   /   \
[1]     [3]

Výsledný diff by měl být maličký – doslova: sem přidej jeden uzel a hotovo.

A právě tady do hry vstupuje Edward Kmett. Vzpomněl jsem si totiž na jeho přednášku, kde mluvil o succinct datových strukturách a wavelet stromech, které mají efektivní operace rank a select. Pak zmínil, že tohle všechno jde použít k práci s JSON dokumenty bez toho, aby je bylo nutné parsovat. Zmínil, že je (za pomoci succint DS) možné vytvořit velice kompaktní index, s jehož pomocí se dají efektivně číst data z JSON stringu bez jeho naparsování do stromové datové struktury. Hlavní ingrediencí tohoto přístupu byly závorky. Ty určovaly, kde daný podstrom začíná a kde končí. V kombinaci s rychlým rank/select pak můžu efektivně najít konec podstromu a celý ho přeskočit.

Závorky jsou hlavní i pro diffování stromů.

Nejdřív musím strom převést na lineární formu, kde závorky udávají začátek a konec každého z podstromů1.

Z dvou příkladů uvedených výše vzniknou následující sekvence tokenů. Číslo je jeden token stejně jako otevírací a uzavírací závorka.

      (5 (2 1 3) 6)
( 4 0 (5 (2 1 3) 6))

Teď když mám dva stringy, můžu spočítat jejich diff přesně tak, jsem to psal minule. Výsledkem bude něco jako:

+(+4+0             -)

A to je všechno. Diff zahrnuje jen nezbytně nutné množství změn: Přidat jeden uzel na začátek s hodnotou 4, levý podstrom obsahuje nulu, pak následuje zbytek stromu beze změny a nakonec je nový strom ukončen. Změny nejsou tak snadno čitelné jako je JSON diff, ale dokážou zahrnout složité transformace a změny struktury.

Ukážu příklad s rotací stromu. Strom před rotací:

        (6)
       /   \
    (2)     [7]
   /   \
[1]     (4)
       /   \
    [3]     [5]

Po rotaci:

      (4)
     /   \
  (2)     (6)
  / \     / \
[1] [3] [5] [7]

Tokenizace, LCS (nejdelší společná subsekvence) a diff:

A:      (  6 (2 1 ( 4 3       5 ) ) 7)
B:      (4   (2 1     3 ) ( 6 5     7))
LCS:    (    (2 1     3       5     7)
diff:   +4-6     -(-4) +)+(+6  -)-)  +)

Je vidět, že diff obsahuje jen uzly 4 a 6. Všechny ostatní zůstaly na svém místě, protože se nezměnilo jejich relativní pořadí.

Po aplikaci diffu už zbývá jen z nové sekvence tokenů postavit nový strom, který obsahuje aplikované změny.

Tohle všechno se dá dělat líně: Sekvenci tokenů nemusím materializovat, můžu ji produkovat líně, stejně tak i patchování může být líné a výsledek je zase líná sekvence ze které líně stavím strom. Pořád je však třeba postavit zcela nový strom.

Dalším logickým krokem je diff interpretovat jako sérii transformací, kterou je možné provést in-place. Ale to až někdy příště, protože zatím nevím, jak to udělat.


Dodatky:

Běžný algoritmus počítající LCS běží v kvadratickém čase a potřebuje kvadratický prostor. Pomocí metody čtyř rusů se dá dostat na čas O(n2 / log(n)). Tedy tak se to aspoň všude píše, ale nikde jsem nenašel veřejně dostupnou publikaci, která by vysvětlila detaily. Jiný algoritmus potřebuje jen lineární místo, další zas běží v čase O((n+r)log(n)), kde r je editační vzdálenost. Speciální případy (např. neopakující se tokeny) se dají vyřešit rychleji. Dále existují techniky aproximace LCS, které jsou rychlejí než přesný výpočet.

Když jsem začal hledat, jak dělat přímo diff stromů bez převedení na sekvenci, v A survey of potential XHTML diff/merge algorithms jsem narazil na XyDiff. Ten běží v běží v průměrném lineárním čase, produkuje přibližné (ale dostatečně dobré) výsledky a je založen na hashování. Každý podstrom je zhashován a to umožňuje rychle hledat odpovídající podstromy v obou vstupech. Tohle by se dalo použít pro snížení počtu tokenů i v mém přístupu.

Kód tady.


  1. Když jde o strom, který má fixní počet potomků, uzavírací závorky nejsou potřeba, protože jsou implicitní.

Flattr this!

Posted in Algo | Leave a comment

diff a komprese

Někdy se zkrátka daří.

Celé dny se nestane nic a pak přijde den, kdy zjistím, jak funguje unixová utilita diff a použiji těchto znalostí, abych komprimoval 3× rychleji a 70× lépe než ctihodný gzip.

Ale začnu popořadě.

Mám dataset, který je tvořen JSON soubory. Nezáleží kolik přesně jich je, jde o něco jako O(hodně) a v katalogu Dellu jsem nenašel stroj, který by měl RAM tak velikou, že by se do ní všechno vešlo. Nejde jen o náhodné JSON soubory, ale log JSONů (stylem jeden JSON na řádek), které jsou si velice podobné. Záznam může mít něco mezi jedním a sto kilobajty s tím, že se od toho předchozího liší jen změnami na několika málo místech: Nějaké vnitřní vlastnosti mohou být aktualizované, do pole může být přidán jeden nebo dva dílčí objekty a jiný může chybět. Podobnost je veliká a většinou přasahuje 99%.

Když je entropie celého datasetu takto malá, kompresní algoritmy by měly zářit.

Bohužel tomu tak není a Gzip v nejlepším případě komprimoval na čtvrtinu. To mi dlouhou dobu přišlo jako neuspokojivý výsledek s přihlédnutím k obrovské repetici v datech. Dlouho jsem to nijak neřešil, protože disky jsou velké a nestojí skoro nic. Ale data narůstala a narůstala, až do té míry, že jsem o zase začal přemýšlet o lepší kompresi.

Napadlo mě, že když jde o stringy, mohl bych uchovávat jen diffy mezi řádky v logu. To na první pohled vypadá jako jednoduchý problém, který se dá sfouknout za půl hodiny. Když jsem o tom ale začal přemýšlet, došlo mi, že to není tak jednoduché. Naivní algoritmy můžou fungovat v mnoha případech, ale dají se velice snadno zmást a pak nenajdou nejmenší možný diff, ale nějakou monstrozitu, která nic moc neušetří.

Potřeboval jsem něco zásadového, potřeboval jsem algoritmus.

Napadala mě parafráze na Kapitána Willarda z Apokalypsy:

Everyone gets everything he wants. I wanted an algorithm. And for my sins, Wikipedia gave me one. Brought it up to me like room service.

Jak jsem velice rychle zjistil, základem unixové utility diff je algoritmus pro nalezení nejdelší společné podsekvence (longest common subsequence/LCS).

Když mám například dvě sekvence/stringy:

A B D E I J K L X
B C D I J L M N

a správně je zarovnám

A B   D E I J K L     X
  B C D   I J   L M N

je vidět, že jejich LCS je

B D I J L

Když mám tuhle sekvenci, můžu snadno spočítat diff. Z prvního řetězce byly odstraněny znaky A E K X a do druhého přidány C M N. Hotovo.

Řešení LCS se dá najít rekurzivně v nepoužitelně pomalém čase nebo přes dynamické programování s pomocí O(n2) extra prostoru. Paměťové nároky situaci komplikují, protože LCS na úrovni znaků dvou stringů, které mají 100kB, by potřebovalo 40GB pomocnou matici. Z toho důvodu utilita diff nebere v potaz jednotlivé znaky, ale celé řádky, kterých je mnohem méně. To dělá problém lépe zvládnutelný.

Já však neměl text, ale jednořádkové JSONy.

Chvíli jsem začal uvažovat, že bych JSONy naparsoval a počítal jejich diffy na úrovni datových struktur. Ale pak jsem, nejspíš pod vlivem ducha Edwarda Kmetta, dostal nápad. Co kdybych string JSONu rozdělil podle čárek a počítal diff na úrovni těchto segmentů. Velice snadno bych tak dostal bloky textu, které mají v JSONu strukturální význam bez toho, aniž bych ho musel parsovat.

Diff mezi těmito dvěma JSONy

{"a":1,"b":2,"c":3}
{"a":1,"b":47,"c":3}

by vypadal

- 7 6       // odstraní 6 znaků počínaje pozicí 7
+ 7 "b":47, // na pozici 7 vloží nový řetězec

To funguje, ale proto, že JSON obsahuje hodně věcí oddělených čárkou a algoritmus běží v podstatě v kvadratickém čase, není to tak rychlé, jak by to jen mohlo být. Pokud JSON obsahuje vnořené objekty, můžu ho nasekat podle uzavírací složené závorky }1. Těch bývá méně a hledání LCS je pak rychlejší.

Ok, teď tedy vím, jak udělat diff JSONů. Ale je to k něčemu dobré?

Místo vysvětlování ukážu malou tabulku:

jsonlog            2926466947 B   originální soubor
jsonlog.gz          772794031 B   komprimováno gzipem
jsonlog.diffs        80895070 B   komprimováno jako série diffů
jsonlog.diffs.gz     12725864 B   diff + gzip

Jak je vidět, tak v tomto případě je samotné převedení na sérii diffů skoro 10× efektivnější než gzip a kombinace s gzipem vede ke 230× kompresi, velice blízko ke skutečné informační hodnotě datasetu.

Kdyby se někdo zajímal, proč si gzip vede tak zle, tak se musím přiznat, že nevím. Moje teorie je, že nevidí, a tedy ani nedokáže využít, podobnosti ve velkém měřítku a komprimuje v malém pracovním okně. Moje doménové znalosti v tomto případě fungují lépe než obecný gzip2. Uvažoval jsem také o xz. Ten komprimuje o něco lépe, ale je strašlivě pomalý.

Samotné diffování je 3× rychlejší než gzip a diffy se na mém stroji (v jednom vlákně) dekomprimují rychlostí 336 MB/s, což je víc než rychlost čtení/zápisu rotujícího disku. Vyplatí se tedy online dekomprimovat data než je číst z disku.

A to podle mě, stojí za tu námahu.


Zdrojový kód je v tomto gistu a tomto repozitáři


Dodatek #1: Ukazuje se, že zstd s maximální úrovní komprese komprimuje stejně dobře jako můj diff+gzip hack a dekomprese je velice rychlá.

Dodatek #2: Po chvíli testování zstd vypadá jako ideální program na obecnou komprimaci. Úrovní komprese je velice blízko xz/lzma, ale pracuje mnohem rychleji a dokáže data dekomprimovat rychlostí přes 1GB za sekundu.


Poznámky:

  1. Pokud na to chcete jít o něco chytřeji a vyhnout se složeným závorkám v řetězcích, stačí udělat dvoustavové parsování: mimo string a ve stringu. Když jsem mimo string, závorky beru v potaz, když jsem ve stringu, přeskakuji je.
  2. Funguje to takhle dobře jen pokud je podobnost mezi sousedícími záznamy veliká. Když by můj dataset byl tvořen zcela náhodnými JSONy, diffová komprese by neušetřila skoro žádné místo, narozdíl do gzipua, který by stále dokázal srazit velokost na něco kolem jedné čtvrtiny.

Flattr this!

Posted in Algo | 8 Comments