Cum se rezolva?Problema este tratata in revista Doctor Dobbs' Journal
(DDJ), o revista apreciata din SUA, cu un continut eterogen (de la
algoritmica la probleme serioase de ingineria programarii). Aceasta
problema a fost abordata chiar in 2 numere din aceasta revista
(aprilie 1998 si aprilie 1999).Sursa Ralucai Sauciuc este izbitor de
asemanatoare cu pseudocodul prezentat intr-unul din articolele din
DDJ, ceea ce ma duce cu gandul ca trebuie sa fi avut respectivul numar
din DDJ. Asadar, sursa oficiala a acesteu runde o voi decreta pe cea a
Ralucai. Asta nu inseamna ca sunt scutit de a explica cum anume
functioneaza algoritmul.De asemenea sursa lui Florin Ghetu pare sa fie
o implementare similara cu algoritmul implementat de Raluca.
Despre algoritm
---------------
Vom porni de la algoritmul forta bruta aplicat acestei probleme si vom
refina treptat-treptat o varianta de complexitate-timp O(m*n) si
complexitate-spatiu O(m).
Algoritmul forta bruta functioneaza in felul urmator: formam toate
sub-matricele continute in matricea data si verificam pentru fiecare
din aceste matrice daca e compusa doar din elemente de 1. Daca da,
verificam daca aria sub-matricei curente este mai mare decat aria
maxima gasita pana acum. Daca da, tinem minte aceasta noua arie drept
aria maxima gasita pana acum.Matricea se formeaza prin alegerea
elementelor extremitati ll (lower-left) si ur (upper-right).
Atentie: prin conventie spunem ca un element ll este lower left de alt
element ur daca llx
O(m*n).
Asadar complexitatea totala a algoritmului forta-bruta este produsul
celor 2 complexitati calculate, adica O(m^3*n^3). Sursa care
implementeaza forta-bruta este Comisia\M3_N3\ver1.cpp.
Foarte ineficient, nu-i asa? Putem imbunatati acest algoritm observand
ca daca o sub-matrice contine un element 0, atunci toate matricele
care contin aceasta sub-matrice vor contine si ele respectivul element
0. Sursa acestei implementari este Comisia\M2_N2\ver2.cpp.
Vom proceda astfel: pentru fiecare varianta posibila a extremitatii ll
ne vom extinde (a se vedea functia GrowOnes) cu extremitatea (x,y)
(extremitatea (x,y) este candidatul pentru (urx,ury) care maximizeaza
aria pentru (llx,lly) curent) doar pe portiunile compacte cu
1. Aceasta extindere se face luand initial y <- lly, si x <- llx,
extinzandu-ne la dreapta pana cand intalnim 0. Cand intalnim 0 tinem
minte valoare coloanei imediat la stanga de coloana curenta in
xMax. Pe urma urcam pe un nivel mai sus si repetam procesul de
extindere, avand grija sa nu ne extindem la dreapta mai mult de
xMax.Sa luam un exemplu:
- avem matricea
c c
o o
l l
0 3
linia 6 0 0 1 1
1 1 1 1
1 1 1 1
1 1 1 0
1 1 0 0
1 1 1 1
linia 0 1 1 1 0
unde elementul cel mai stanga jos este elementul (0,0), iar
elementul cel mai de dreapta-sus este (6,3).
- avem llx=0, lly=0. A se consulta in continuare funtia GrowOnes.
- initial avem xMax=3. Pornim cu y=0 si x=0 si ne extindem spre
dreapta. Mergem pana la ultimul element 1 de pe aceasta
linie. Acest element este y=0, x=2. Facem xMax=2 si, verificam daca
aria sub-matricei (0,0)->(0,2) este mai mare decat aria maxima
gasita pana acum (adica 0) si facem maxArea=3.
- urcam cu un nivel (il incrementam pe y --> y=1). Ne extindem din
nou spre dreapta, tinand minte ca nu ne putem extinde mai mult de
coloana 2. Avand doar 1 pe aceasta linie, xMax ramane
neschimbat. Verificam daca aria sub-matricei (0,0)->(1,2) este mai
mare decat aria maxima gasita pana acum --> facem maxArea=6.
- urcam cu un nivel (il incrementam pe y --> y=2). Ne extindem din
nou spre dreapta, tinand minte ca nu ne putem extinde mai mult de
coloana 2. Intalnim 0 de pe coloana x=2. Il facem pe xMax <-
1. Verificam daca aria sub-matricei (0,0)->(2,1) este mai mare
decat aria maxima gasita pana acum, adica 6. Nu este mai mare, asa
ca nu modificam maxArea.
- si tot asa pana sus, functia GrowOnes asigurandu-ne gasirea
maximului ariei sub-matricei cu una din extremitati ll. Acest
algoritm are o complexitate de O(m^2*n^2).
In continuare mai putem sa aducem imbunatatiri, observand ca facem mai
multe extinderi la dreapta asupra acelorasi elemente in cadrul
aceleiasi executii. Atunci putem sa tinem evidenta intr-un vector (cei
din DDJ ii spun cache) cu m elemente, cate elemente consecutive de 1
sunt de la coloana curenta (llx) pana in marginea dreapta a
matricei. Ceea ce aduce un spor de eficienta este faptul ca se poate
calcula foarte usor (in O(m)) cache-ul coloanei imediat la dreapta din
cache-ul coloanei curente.Sursa acestui algoritm este
Comisia\M2_N2\ver3.cpp.
Sa vedem pe cazul matricei prezentate mai sus care este evolutia
cache-ului. Cache-ul este initializat cu 0. Cache-ul este reprezentat
pe coloane. Coloana cea mai din dreapta (doar cu elemente 0)
reprezinta starea lui initiala. A doua coloana din dreapta reprezinta
starea cache-ului pentru llx=3, a treia coloana starea cache-ului
pentru llx=2, etc.
0 0 2 1 0
4 3 2 1 0
4 3 2 1 0
3 2 1 0 0
2 1 0 0 0
4 3 2 1 0
3 2 1 0 0
Vom adapta functia GrowOnes din versiune precedenta a algoritmului,
astfel incat sa foloseasca valorile acestui cache in loc sa ne
extindem la dreapta pas cu pas pentru fiecare y. Acest algoritm are o
complexitate-timp de O(m^2*n), platindu-se costul maririi
complexitatii-spatiu cu O(m) (datorita introducerii cache-ului). O
observatie demna de mentionat este ca daca avem m>n putem sa mai
reducem timpul de executie aplicand algoritmul asupra transpusei
matricei, aria maximala fiind invariant la transpunerea matricei.
Ce imbunatatiri mai putem aduce pentru a reduce complexitatea-timp la
O(m*n)? Observam ca ultimul algoritm verifica dreptunghiuri care sunt
continute in totalitate de alte dreptunghiuri verificate. Spre exemplu
algoritmul precedent verifica sub-matricele (0,0) -> (0,2) si (0,0) ->
(1,2), cea de-a doua continand-o pe prima. Scopul nostru este sa
examinam doar dreptunghiurile care nu sunt continute in alte
dreptunghiuri. Sa luam un exemplu, in care avem cache-ul pentru
coloana curenta cu urmatoarele valori:
1 2 3 4
linia 9 0
4 . . . .
4 . . . .
3 . . .
2 . .
4 . . . .
3 . . .
1 .
1 .
linia 0 0
Coloana din stanga reprezinta cache-ul. Pentru a intui mai clar se
poate recrea matricea de 0 si 1 corespunzatoare cache-ului descris.
Coloanele din dreapta cu linii verticale (care sunt notate deasupra cu
1, 2, 3, 4) reprezinta cat de mult se intind pe verticala
dreptunghiurile cu 1 care nu sunt continute in alte dreptunghiuri cu
1, de latime 1, 2, 3, respectiv 4, care se intind pe orizontala de la
coloana corespunzatoare cache-ului in dreapta.
Evident coloanele cu liniile verticale se deduc din valorile
cache-ului, astfel:
- pornim cu elementul din cache corespunzator liniei 0 din matrice
si "urcam" in cache, pas cu pas, spre elementul corespunzator
liniei 9.
- daca elementul curent din cache este mai mare decat anteriorul
atunci inseamna ca trebuie sa inceapa una sau mai multe linii
verticala corespunzatoare unor latimi cel mult egale cu elementul
curent
- daca elementul curent din cache este mai mic decat anteriorul
atunci inseamna ca s-au terminat una sau mai multe linii
verticale corespunzatoare unor latimi cel putin egale cu
elementul curent
- daca elementul curent din cache este egal cu anteriorul nu incepe
si nu se termina nici o linie verticala.
Se observa ca liniile verticale au o cvasi-proprietate a parantezelor
bine formate: o linie verticala care incepe mai devreme se termina
intotdeauna mai tarziu sau in acelasi moment cu o linie care a inceput
mai tarziu; in cazul in care 2 linii verticale pornesc la acelasi timp
atunci cea corespunzatoare unei latimi mai mari se termina prima.
Aceasta proprietate ne duce cu gandul sa folosim o stiva in care:
- facem PUSH atunci cand elementul curent din cache este mai mare
decat cel anterior, cu perechea de valori
(valoarea_elementului_anterior; pozitia_lui_in_cache)
- facem POP-uri in cazul in care elementul curent este mai mic
decat cel anterior, pana cand extragem o pereche care are
valoare_element mai mica decat elementul curent din c. Pentru
fiecare POP fac verificarea daca dreptunghiul de latime
valoare_element si inaltime
y_curent-pozitie_in_cache_scoasa_din_stiva are aria mai mare
decat aria maxima gasita la pasul curent.
Sursa acestui algoritm este Comisia\M_N\sauciuc.cpp. Sa analizam
complexitatea-timp a algoritmului: se executa de m ori functia care
face actualizarea cache-ului (care functie are complexitate Theta(m)),
iar bucla for(j=0; j<=n; j++) are complexitate O(n), pentru ca se fac
cel mult n+1 PUSH-uri, iar numarul de POP-uri este egal cu cel de
PUSH-uri. Complexitatea spatiu este tot O(n), chiar daca s-a mai
introdus o stiva (care are dimensiune n+1).
Evident aici ia sfarsit questul nostru de reducere a complexitatii,
din moment ce este nevoie de exact m*n operatii de citire a
elementelor din matrice, iar daca nu sunt citite toate este posibil sa
nu se gasesca solutia optima. Pe deasupra aceasta varianta finala de
algoritm este posibil sa o implementam fara a tine minte in memorie
toate matricea (in celelalte versiuni ale algoritmului trebuia sa
retinem intreaga matrice in memorie). Mai exact nu trebuie sa retinem
nici un element din matrice: citim element cu element din linia
curenta si actualizam cache-ul. Practic, reducand complexitatea-timp
am redus-o si pe cea spatiu.
Se observa ca in program (Comisia\M_N\sauciuc.cpp) se citesc
elementele linie cu linie, deci cache-ul se calculeaza pe matricea
transpusa fata de cum a fost prezentat algoritmul mai sus.
Un alt tip de algoritm
----------------------
Spre norocul tuturor, mai exista un algoritm care rezolva problema in
O(m*n). Acesta este dupa parerea mea mai intuitiv decat cel anterior,
dar are o implementare un pic mai costisitoare. Programul care
implementeaza algoritmul pe care il voi descrie mai departe este "M_N
ver 2\algo2.cpp".
Care sunt pasii acestui algoritm:
- se citeste matricea cate o linie la un pas
- se calculeaza in vectorul C, cache-ul despre care am discutat in
algoritmul anterior (adica c[i]=numarul consecutiv de 1 care sunt
incepand de la linia curenta, deasupra).
- avem acum un vector C cu n elemente. Cum anume putem sa
determinam la pasul curent, in timp O(n) care este aria maxima
pentru pasul curent? Pentru asta trebuie sa gasim lungimea unui
interval de elemente din C, (C[i1..i2] - interval solutie) care
sa maximizeze aria, adica valoarea:
(i2-i1+1) * min(C[i])
i1<=i<=i2
Fie M elementul maxim din C. Sa observam ca daca M apartine
intervalului solutie, avem 2 variante: fie intervalul solutie este
compus doar din elementul M (atunci aria maxima este 1*M=M), fie
intervalul contine si un element vecin (poate chiar mai multe) al lui
M (care este mai mic sau cel mult egal cu M), in care caz aria maxima
se obtine dubland valoarea elementului vecin lui M. Spunem in cazul
asta ca elementul vecin lui M contribuie de 2 ori cu valoarea lui, in
calculul ariei.
Definim vectorul length de n elemente, unde length[i] inseamna cati
vecini luam in calcul cand luam in consideratie elementul i din
vectorul C pentru calcul ariei.
Sper ca se vor clarifica lucrurile imediat. Iata cum determinam in
O(n) aria maxima:
- bagam elementele lui C bagate intr-o lista L
- initializam elementele vectorului length cu 1
- atata timp cat L are elemente:
P1. gasim maximul M din lista L. (pentru a obtine complexitate
O(n) pe toata aceasta bucla, se recomanda a se face o
sortare initiala a tuturor elementelor lui C cu un
counting-sort care are complexitate O(n), a.i. sa alegem
pe urma la fiecare pas maximul din L in O(1) - adica al
k-lea maxim, unde k=nr pas)
P2. verificam daca length[i]*M (unde length[i] este elementul
lui length asociat lui M) este mai mare decat cea mai buna
arie gasita pana acum, in care caz o retinem
P3. incrementam length[vecinul_cel_mai_mare_al_lui_M_din_L] cu
length[i]
P4. stergem M din L.
Sa luam un exemplu. Fie C=(4,2,3,6,5,3).
- atunci la primul pas avem:
length=(1,1,1,1,1,1)
L=(4,2,3,6,5,3)
Gasim maximul din L --> M=6. Verificam daca 1*6>maxArea -->
maxArea=6. Gasim vecinul cel mai mare al lui 6 din L (acesta este
5) caruia ii incrementam length cu length-ul lui 6.
- La al doilea pas avem:
length=(1,1,1,1,2,1)
L=(4,2,3, ,5,3)
Gasim maximul din L --> M=5. Verificam daca 2*5>maxArea -->
maxArea=10. Avem ambii vecini din L ai lui 5 egali, il alegem pe
oricare, spre ex. 3-ul din dreapta.
- La al treilea pas avem:
length=(1,1,1,1,2,3)
L=(4,2,3, , ,3)
Gasim maximul din L --> M=4. Verificam daca 1*4>maxArea -->
maxArea=10. 4 are un singur vecin...
- La al 4-lea pas avem:
length=(1,2,1,1,2,3)
L=( ,2,3, , ,3)
Gasim maximul din L --> M=3 (din dreapta). Verificam daca
3*3>maxArea --> maxArea=10. 3 are un singur vecin...
- Pasul 5:
length=(1,2,4,1,2,3)
L=( ,2,3, , , )
Gasim maximul din L --> M=3. Verificam daca 4*3>maxArea -->
maxArea=12. 3 are un singur vecin...
- Pasul 6:
length=(1,6,4,1,2,3)
L=( ,2, , , , )
Gasim maximul din L --> M=2. Verificam daca 6*2>maxArea -->
maxArea=12. L devine vid, deci se termina executia algoritmului.
Observam ca faptul ca ajungem la optim este garantat de ordinea in
care se verifica intervalele candidat.
Obs: In implementare nu respect pasul P1 al algoritmului, anume nu
gasesc maximul din L (trebuia implementat ceva mai mult), ci gasesc un
maxim local, adica un element care are ambii vecini mai mici sau cel
mult egali cu el. Daca aveti timp si chef, va invit sa va uitati la
program si sa incercati sa generati un caz care sa imi faca programul
sa aiba timpul de executie mai mare de O(m*n). Eu unul nu am fost in
stare sa generez.
Se observa ca in cazul algoritmului cu stiva se pot determina si care
sunt coordonatele dreptunghiului solutie imediat, pe cand in cazul
algoritmului cu lista dublu-inlantuita aceasta inseamna a cauta in
vectorul length, intervalul de lungime length[i] care sa contina
pozitia i si care sa contina toate valorile 1..length[i] o singura
data.
Concluzie: din culpa am dat niste solutii kilometrice. Daca am sa mai
propun probleme promit solemn ca le voi rezolva in timp util si cat
mai simplu si cat mai putin cu putinta.
Referinte:
1. David Vandervoorde - The Maximal Rectangle Problem, DDJ, April 1998
2. Udi Manber - Designing Algorithms Incrementally, DDJ, April
1999
|
|
