Często problem złożoności obliczeniowej można odłożyć na później. Bo działa. Ale kiedy trafi się, że nasz algorytm będzie musiał przerobić większe ilości danych, to nagle to czy działa liniowo, kwadratowo, a może i sześciennie ma ogromne znaczenie i jest różnicą między otrzymaniem wyniku za 10 minut lub za 10 dni.
Gdzie spotykamy się z dużymi danymi? Na niektórych konkursach programistycznych, a przede wszystkim - w prawdziwym życiu (Big Data).
Zacząłem w tym semestrze przedmiot Wstęp do biologii obliczeniowej. Brzmiał ciekawie, więc się zapisałem. Lekki on dla mnie nie jest, bo przez ostatnie 5 lat z biologią styczności nie miałem. W dodatku mam kolizję wykładu z WBO z ćwiczeniami z Zaawansowanego programowania funkcyjnego. Skończyło się tym, że czytam źródłową książkę, na podstawie której opracowano przedmiot.
Na pierwszych zajęciach został przedstawiony problem sekwencjonowania DNA. W skrócie: w laboratorium, żeby badać DNA to trzeba go mieć dużo, więc się je klonuje przy użyciu bakterii. Następnie się używa pewnego enzymu, który kroi DNA na kawałki. Ale może pokroić je w różnych miejscach, więc trzeba rozwiązać problem odbudowania oryginalnej sekwencji z wielu podsekwencji.
I tu wchodzą algorytmy. Można tu znaleźć problem znajdywania ścieżki Hamiltona (przez wszystkie wierzchołki raz) w grafie zachodzących na siebie kawałków DNA. Ten problem jest NP, ale okazuje się, że lekko modyfikując graf, można zejść do problemu znajdywania ścieżki Eulera (przez każdą krawędź raz), który ma sensowny algorytm.
Ale ja piszę ten post, żeby opowiedzieć o moim zadaniu domowym.
Mieliśmy zaimplementować algorytm Pythonie, który:
- Wczyta wybrane sekwencje DNA drożdży.
- Dla każdej sekwencji DNA znajdzie unikalny k-probe, który odpowiada jakiejś podsekwencji.
- Sprawdzi jakie jest najmniejsze k, dla którego zbiór unikalnych k-probe istnieje.
Sekwencji drożdży było jakieś 8,8tyś, z czego najkrótsza sekwencja miała 702 nukleotydy.
Żeby wygenerować wszystkich kandydatów na probe użyłem takiej funkcji
def spectrum(s, k):
for begin in range(0, len(s) - k + 1):
yield s[begin:begin+k]
Następnie przygotowałem implementację funkcji sprawdzającej czy dla danego k da się znaleźć ten zbiór existsUnique(dna, k). Zaraz ją omówię.
Następnie napisałem funkcję, która dostaje długość najkrótszej sekwencji, czyli maksymalne k i przeprowadza wyszukiwanie binarne, żeby znaleźć najmniejsze k dla którego existsUnique zwróci True.
Dobra, więc sedno tego posta, to algorytm existsUnique.
Początkowo próbowałem to zrobić na zbiorach
spectrs = [set(spectrum(elem.seq, k)) for elem in dna]
for idx in range(0, len(spectrs)):
temp = spectrs[idx]
for odx in range(idx + 1, len(spectrs)):
spectrs[idx] = spectrs[idx] - spectrs[odx]
spectrs[odx] = spectrs[odx] - temp
return allNonEmpty(spectrs)
Ale jest to O(len(spectrs)^2 * maxlen(spectrum)) i jak zostawiłem komputer żeby to robił przez noc, to 8,5h minęło, a on nie sprawdził nawet jednego k = 351.
Potem spróbowałem podejścia od innej strony (combinations zwraca po kolei wszystkie k-słowa z alfabetu DNA)
spectrs = [set(spectrum(elem.seq, k)) for elem in dna]
P = [None for s in spectrs]
for probe in combinations(k):
counter = None
for idx,s in enumerate(spectrs):
if probe in s:
if counter != None:
counter = None
break
counter = idx
if counter != None:
P[counter] = probe
if allNonEmpty(P):
print(P)
return True
return False
To rozwiązanie również jest strasznie wolne ~ O(4^k * len(spectrs) * log(maxlen(spectrum)).
Więc zacząłem się zastanawiać nad strukturami danych jakich używałem. Spróbowałem użyć struktury Trie
class TrieNode(object):
def __init__(self, char):
self.char = char
self.children = [None, None, None, None]
self.seqs = -1
def isLeaf(self):
return (self.children[0] == None and self.children[1] == None
and self.children[2] == None and self.children[3] == None)
# num(c) A => 0, T => 1, G => 2, C => 3
class Trie(object):
def __init__(self):
self.root = None
def insert(self, word, idx):
node = self.root
for char in word:
if node:
if node.children[num(char)]:
node = node.children[num(char)]
else:
newnode = TrieNode(char)
node.children[num(char)] = newnode
node = newnode
else:
self.root = TrieNode(char)
node = self.root
# node = leaf of the inserted path
# add the id of DNA sequence to the set
# of sequences containing this subsequence
if node.seqs == -1:
node.seqs = idx
else:
node.seqs = -2
# walk over every leaf and check its
# set of sequences. if it has only one
# seq then it can be its probe
# add that seq will appear in the returned set
def singletonSet(self):
self.seqs = set()
def walk(node):
if node.isLeaf() and node.seqs > 0:
self.seqs |= node.seqs
else:
for c in node.children:
if c:
walk(c)
if self.root:
walk(self.root)
return self.seqs
# for a given length `k` of probes
# we want to try to find the set P
def existsUnique(dna, k):
print("Trying k = ", k)
t = Trie()
for idx, seq in enumerate(dna):
for s in spectrum(seq.seq, k):
t.insert(str(s), idx)
# if for every dna seq it belongs to
# the singletonSet of trie then
# there exists a desired set P
return len(t.singletonSet()) == len(dna)
Algorytm jest liniowy od podsekwencji, razy logarytm z k. Ma jedną wielką wadę - pożera niesamowicie pamięć. Tak bardzo, że nie byłem w stanie go uruchomić dla k > 15.
Zacząłem myśleć - może by tak wrócić do najprostszych rozwiązań?
Weźmy słownik, liniowo się przejdźmy po spektrach i dodajmy każde do słownika. Jak wcześniej takiego nie było to dodam indeks sekwencji w której on występuje, a jak był to ustawię liczbę ujemną. Na koniec zbiorę indeksy do setu i sprawdzę, czy jest ich tyle co sekwencji.
def existsUnique(dna, k):
print("Trying k =", k)
t = {}
for idx, seq in enumerate(dna):
for s in spectrum(seq.seq, k):
if not str(s) in t:
t[str(s)] = idx
# if t[str(s)] == idx it means
# we encountered a duplicate
# subsequence in the same seq
elif t[str(s)] != idx:
t[str(s)] = -2
# probable is the set of sequences
# that have a probe
probable = set([x for x in t.values() if x >= 0])
# if all sequences have a probe then
# a set P exists
return len(probable) == len(dna)
I tak oto cały program drastycznie przyśpieszył i zakończył się po 2,5 minutach.
I to koniec tej pięknej historii szukania optymalnego algorytmu.