Transcript 1. Inledning

1. Inledning
Inom fylogenin, dvs läran om arters utveckling och deras inbördes släktskap,
har man redan länge studerat organismer och bildat släktträd (kladogram)
mellan dessa utgående från arternas yttre form och struktur. I och med
molekylärbiologins
framsteg
har
nya
tekniker
gjort
analysering
och
dokumentering av organismers genetiska material (DNA) möjliga. Genom att
tillämpa DNA-data inom fylogenin har det uppstått helt nya möjligheter att
klassificera organismer och att konstruera kladogram. Eftersom genetiska det
materialet innehåller data som beskriver organismers utseende och beteende
är det naturligt att dessa data kan användas vid analys och uppbyggnad av
kladogram.
DNA-sekvenserna som används i fylogenetiska analyser är ofta långa, varför
det krävs datorprogram för att kunna generera pålitliga resultat. POY
(Phylogenetic Analysis of DNA and other Data using Dynamic Homology), se t.ex.
[Whee06 s. 29 – 35], är ett datorprogram som läser in data i form av DNAsträngar och genom analys och manipulering av dessa producerar kladogram
som resultat. POY grundar sig på ett koncept som kallas för dynamisk
homologi
vilket
innebär
att
två
tekniker
(strängoptimering
och
kladogramsökning) knyts samman. Strängoptimering innebär att DNAsträngar för par av organismer jämförs med hjälp av en teknik som kallas
inpassning. I kladogramsökning utnyttjas resultatet som producerats av
strängoptimeringarna (i form av ”distanskostnader” mellan DNA-strängarna)
och
olika
heuristiska
sökmetoder
används
för
att
åstadkomma
verklighetstrogna kladogram.
Processen att skapa ett kladogram kan, beroende på indata, dröja dagar,
veckor eller t.o.m. månader att utföra, varför en uppsnabbning av programmet
spelar en central roll. I denna avhandling kommer olika kodoptimeringar att
göras för att minska programmets körtid.
1
1.1. Avhandlingens ändamål
Ett av syftena med denna avhandling är att försöka ge läsaren en översikt av
hur det kladistiska programmet POY fungerar. Genom att ta upp koncept som
strängoptimering och kladogramsökning ska läsaren förhoppningsvis kunna
bilda sig en uppfattning om hur kladogramkonstruktion går till.
Det viktigaste målet med avhandlingen är att försöka visa hur tillämpning av
olika optimeringsmetoder kan snabba upp exekveringen av ett program; i
detta fall POY. Genom att bl.a. presentera optimeringsmetoder som
eliminering av hopp och upprullning av slingor ska läsaren ges en bild av hur
mycket snabbare ett redan tämligen optimalt program kan göras.
1.2. Strukturering och uppläggning
Denna avhandling inleds i kapitel två med att ge förklaringar till begrepp som
fylogeni och kladistik. Läsaren kommer i detta kapitel att ges information om
vad ett kladogram är och hur organismer ordnas i ett sådant.
I kapitel tre berörs enskilt konceptet med sekvensinpassning. Olika metoder
tas upp för hur man kan ordna sekvenser av DNA så att de sinsemellan kan
jämföras och betraktas på ett objektivt sätt. En viktig metod för att göra
sekvensinpassningar är Needleman-Wunsch algoritmen som kommer att
behandlas närmare. Att sekvensinpassning har placerats i ett skilt kapitel
beror på att det är en av de mest centrala delarna i teorin bakom POY.
I kapitel fyra behandlas POY som en helhet. De två viktigaste delarna av
dynamisk
homologi
(som
är
huvudidén
bakom
POY)
består
av
strängoptimeringsmetoder och kladogramsökningsmetoder. Främst direkt
optimering, som möjligen är den viktigaste strängoptimeringsmetoden,
kommer att tas upp.
2
I kapitel fem förklaras först ytligt hur moderna processorer fungerar. Efter det
presenteras de optimeringsalgoritmer som under arbetets gång utnyttjats och
tillämpats på delar av POY.
I kapitel sex presenteras utdata och optimeringsresultat för olika körningar av
POY.
Även
hårdvarumätningar
som
gjorts
med
hjälp
av
programmeringsgränssnittet PAPI kommer att läggas fram.
3
2. Evolutionsbiologi
Grovt kan en DNA-molekyls sekvens sägas bestå av nukleotider vilka
innehåller fyra typer av nukleotidbaser: Adenin (A), Guanin (G), Cytosin (C)
samt Tymin (T) [Hand07]. Dessa baser sammanfogar, likt figur 2.1, två långa
kedjor av vilka mindre grupper kan delas in i så kallade gener som bestämmer
organismers utseende och beteende.
Figur 2.1. Förenklad DNA-sekvens.
När individer förökar sig kopieras generna till avkomman som för arvet
vidare. Men kopieringen sker inte alltid rätt och i sällsynta fall uppstår en
mutation, vilket innebär att en bas eller en sammanhängande sträng av baser
raderas, byter plats eller införs i DNA-sekvensen. Oftast ger mutationerna
upphov till icke önskvärda egenskaper, men i sällsynta fall kan hos
avkomman uppstå en ny egenskap som är till fördel i den omgivande miljön.
Om t.ex. en flugas DNA muterar så att bakdelen blir mörkare eller så att
flugan blir mera skygg, kan överlevnadssannolikheten öka. När en mutation
äger rum till följd av kopieringsfel under celldelningen förgrenar sig
avkomman från anfadern och en ny undergrupp bildas.
2.1. Fylogeni
Inom
fylogenin
försöker
man
återspegla
den
evolutionära
utvecklingssekvensen för ett antal valda organismer genom att hypotetiskt
gruppera dessa enligt ursprung och bilda överskådliga släktskapsträd
[WSBF91, s.1 – 3]. Detta utförs genom att studera släktskap baserade på nu
levande organismers DNA-sekvenser eller genom att studera organismers
yttre form och struktur. Ett problem vid analys av DNA-sekvenser är att
4
nukleotidmaterialet för döda organismer relativt snabbt bryts ned. På grund av
detta fenomen baseras fylogenetiska studier av DNA främst på levande
organismer.
2.2. Kladistik och kladogram
Kladistik, som också kallas fylogenetisk systematik, är en systematisk
vetenskap som används vid jämförelse av genetiska data i syfte att gruppera
arter utgående från deras antagna evolutionära anfäder. Kladistiska metoder
tillämpas oftast inom fylogenin, men används också inom lingvistiken
[NoCa77]. I denna avhandling nämns kladistik endast i kombination med
fylogenin.
Inom fylogenin försöker man åskådliggöra organismers släktskap sinsemellan
genom att konstruera kladogram som även ibland benämns fylogenetiska träd.
Ett kladogram motsvarar ett fullständigt binärt rotat träd T för vilket, enligt
definition [PrEn97 s.5 – 6], dessa krav bör gälla:
- Det finns inga cykler i T;
- Varje intern nod har en förfader och två avkommor, jämfört med en
lövnod som endast har en förfader, se figur 2.2.1.
- Det finns endast en nod kallad rot, vilken är ett specialfall av en intern
nod i trädet pga. att den inte har någon förfader.
Figur 2.2.1. Nodhierarkin i ett binärt träd.
Varför just binära träd oftast tillämpas inom evolutionsbiologin beror på att
de beskriver mutationskonceptet på ett naturligt sätt. Om avkomman för en
organism har muterat till följd av felaktig kopiering av DNA kommer
5
avkomma och förfader att förgrenas till skilda underträd, vilket i ett binärt
träd kan representeras som en förgrening av en nod. Binära träd är dessutom
enklare att representera i en dator än vad mera generella träd är och enligt
Press m.fl. är det ofta lättare att påvisa olika teorem med hjälp av binära träd
[PTVF07 s. 870].
Lövnoderna i ett kladogram kallas för observerade taxa (Operational
Taxonomic Unit, OTU, se [Whee06 s. 16]) pga. att de innehåller data som
representerar vissa grupper av individer eller organismer (t.ex. människa eller
ko). Lövnoderna är de enda noder i kladogrammet som bär på objektiv
information. Denna information kan bestå av DNA-sekvenser som tagits fram
i laboratoriemiljöer eller av data beskrivande organismernas utseende och
struktur. En fördel med att endast organismerna presenteras som lövnoder är
att det på förhand går att räkna ut hur många interna noder (Hypothetic
Taxonomical Units, HTU, se[Whee06 s. 16]) ett träd kommer att ha [PTVF07
s. 871]. Om ett binärt rotat träd har n OTU:n (lövnoder), har det exakt n - 1
HTU:n (interna noder).
De interna noderna, HTU:na, representerar hypotetiska anfäder i evolutionen
och är endast uppskattningar. Då kladogram konstrueras fungerar dessa som
ett mellanresultat vid sammanknytningar av OTU:n. Genom kombination av
OTU:n med HTU:n och HTU:n med HTU:n uppstår också nya hypotetiska
anfäder.
6
Figur 2.2.1. Exempel på ett kladogram.
Figur 2.2.1 är ett exempel på hur ett kladogram kan se ut. Organismer som är
nära besläktade bör också vara nära varandra i kladogrammet. Även om de är
nära ordnade behöver det nödvändigtvis inte betyda att arterna är nära
besläktade eftersom alla de organismer som testas kan vara väldigt olika. Ett
kladogram beskriver inte en absolut grad av likhet individerna emellan, utan
ställer dem endast i relation till varandra. Om man för flugan, myggan och
människan jämför DNA-sekvenserna för gener av samma typ, är myggan och
flugan närmast besläktade, men ett kladogram ger ingen information om i
vilket skede den gemensamma förfadern till dessa muterade. Förmodligen kan
man heller inte veta hur en exakt förfader till myggan och flugan såg ut. Man
kan inte heller avgöra om t.ex. myggan är en direkt förfader till flugan. Om så
vore fallet kommer de båda ändå att betraktas som syskon i ett kladogram,
vilket kan ställa till problem i en kladistisk analys. Detta beror på att
observerade taxa inte tillåts bli placerade i de inre noderna.
När kladogram konstrueras utgår man från en mängd, ofta oordnade,
organismer vars DNA-sekvenser har studerats i laboratoriemiljöer och
information för dessa erhållits och dokumenterats på ett för datorer
7
tillgängligt sätt (inom datavetenskapen representeras oftast DNA-sekvenser
som strängar av tecken ur ett givet alfabet, där ett tecken motsvarar
begynnelsebokstaven för namnet på en nukleotidbas). Dessa observerade taxa
kan analyseras med hjälp av olika applikationer dedikerade för ändamålet.
Även med hjälp av modern datorkraft kan uppbyggnaden av kladogram vara
en krävande uppgift. Att analysera och ordna hundratals organismer med
DNA-sekvenser på tusentals element kan ta veckor eller t.o.m. månader att
utföra, beroende på hur noggranna sökningarna är.
8
3. Inpassning av DNA-sekvenser
Sekvensinpassning (sequence alignment) eller inpassning kallas den process
där två (eller flera) DNA-sekvenser ställs upp mot varandra, omordnas och
jämförs. Genom att inpassa organismers DNA-strängar kan man lättare avgöra
graden av likhet för dessa. På basis av resultaten av jämförelserna kan
uppbyggnad av kladogram ske med hjälp av olika tekniker.
Ett problem vid jämförelse av DNA-strängar är att utföra själva inpassningen
så att resultatet kan betraktas som pålitligt. En direkt jämförelse av två
motsvarande sekvenser av DNA, bas för bas, vore naiv pga. att DNA vid
mutation inte kanske enbart byter nukleotidbas (substitution), utan möjligtvis
utsätts för en radering (deletion) eller en insättning (insertion) av en eller
flera nukleotidbaser vilket medför att hela strängen kommer att skifta,
position för position [Tola01].
S1: AGGGCTCTTGCTGTCTCTCTCT
||||||||||||||||||||||
S2: AGGGCTCTTGCTGTCTCTCTCT
Exempel 3.1.
I exempel 3.1 är två strängar, S1 samt S2, inpassade mot varandra. I detta fall
är
sekvenserna
identiska
där
”|”
betecknar
likhet
mellan
par
av
nukleotidbaser.
S1: AAAACTCTTGTTGTCTCTCTCT
|
|||||| |||||||||||
S2: AGGGCTCTTGCTGTCTCTCTCT
Exempel 3.2.
I exempel 3.2 har två lokaliserade mutationer uppstått. I S1 har GGG
förändrats till AAA samt ett C har bytts ut mot ett T. I detta sällsynta fall
lyckas en naiv inpassning av sekvenserna eftersom positionerna för de övriga
baserna inte har ändrat.
9
S1: AGGCTCTTGCTGTCTCTCTCT
|||
|
S2: AGGGCTCTTGCTGTCTCTCTCT
Exempel 3.3.
I exempel 3.3 har en radering uppstått i S1 (eller möjligtvis en insättning i
S2). Detta innebär att strängarna som före raderingen matchade varandra
perfekt, nu har blivit nästan helt olika. En bas-mot-bas inpassning i denna
situation skulle ge ett förkastligt resultat, även om strängarna nästan är
identiska. Detta problem kan lösas genom att införa en lucka ”-” i strängen så
att de övriga relativa positionerna bevaras.
S1: A-GGCTCTTGCTGTCTCTCTCT
| ||||||||||||||||||||
S2: AGGGCTCTTGCTGTCTCTCTCT
Exempel 3.4.
I exempel 3.4 har S1 från exempel 3.3 justerats genom insättning av en lucka
i position 2. På detta sätt bevaras den fortsatta ordningen och strängarna
skiljer sig endast vid ett par av baser. Ett icke-trivialt problem vid insättning
av luckor uppstår om en radering i en delsekvens av lika baser har skett. I
exempel 3.4 kunde luckan lika väl införas i position tre eller fyra i S1 och en
sekvensinpassning skulle fortfarande producera samma resultat. Detta betyder
att det existerar tre olika strängar S1’ S1’’ och S1’’’ vilka alla borde tas i
betraktande vid en fortsatt analys.
För att kunna avgöra i hur hög grad två sekvenser överrensstämmer används
olika typer av sekvensinpassningsalgoritmer. Sådana används för att hitta den
bästa sekvensinpassningen mellan två (eller flera) strängar. Ofta används
olika typer av poängräknings- eller kostnadssystem för att kunna avgöra
graden av likhet mellan strängarna vilka ofta är representerade som sekvenser
av tecken eller bitmängder (bit sets). Needleman-Wunsch, se [NeWu70], är ett
exempel på en sådan algoritm och var den första i sitt slag att utnyttja
konceptet med dynamisk programmering (dynamic programming) för att lösa
10
problemet med bästa möjliga sekvensinpassning. Algoritmen beräknar den
optimala inpassningen mellan två strängar genom att bilda en stig mellan övre
vänstra hörnet och nedre högra hörnet en poängmatris F, vars bredd M + 1
och höjd N + 1 motsvarar längderna av strängarna S1 (M) och S2 (N). Initialt
är F(0,0) = F(0, j) = F(i, 0) = 0, där i = 0..M och j = 0..N.
Vidare gäller:
F(i − 1, j − 1) + σ (S1i ,S2 j ) (fall 1)

(fall 2)
F (i, j) = max F(i − 1, j ) − d
F(i, j − 1) − d
(fall 3)

där
d
är
kostnaden
för
en
lucka
i
S1
eller
S2
och
σ
är
en
poängomvandlingsfunktion för inpassning av nukleotidbaserna i positionerna
i respektive j. Fall 1 motsvarar poängtalet efter en inpassning av tecknen
(nukleotidbaserna) S1i och S2 j . Fall 2 och 3 motsvarar poängtalen efter en
horisontell respektive vertikal luckinfogning. Ofta används en pekarmatris P
(också kallad bakåtspårningsmatris) för att hålla reda på den optimala stigen
till varje enskild cell:
P ( i, j ) = P ( i, j ) ∪ diagonal , om fall 1 gäller
P ( i, j ) = P ( i, j ) ∪ vänster , om fall 2 gäller
P ( i, j ) = P ( i, j ) ∪ upp , om fall 3 gäller
Om flera fall gäller samtidigt kommer P(i, j) att peka tillbaka på flera celler
än en. Matrisutfyllnaden avslutas när F (M, N) fyllts. Poängtalet i denna cell
är det optimala resultatet för sekvensinpassningen mellan S1 och S2,
användande poängomvandlingsfunktionen σ. För att reda ut den slutgiltiga
sekvensinpassningen används P vid bakåtspårning. Eftersom varje cell i P
lagrar information om från vilket håll den optimala stigen har kommit, blir
körtiden linjär jämfört med själva poängmatrisfyllnadsproceduren som tar
O(M*N) tid att utföra.
11
Figur 3.1. Fylld poängmatris enligt Needleman-Wunsch, med en utav de bästa
stigarna markerad.
I figur 3.1 har sekvenserna GAATTCA och GGATCG blivit inpassade med
stöd av en uniform poängomvandlingsfunktion, vilket innebär att poängen för
varje omvandling har varit lika stor (i detta fall 1). Bästa poängen för
inpassningen blir i detta fall 4 och pilarna representerar en ur mängden
optimala stigar som i linjäriserad form efter bakåtspårning kan skrivas som:
G-AATTCAGGA-T-C-G
Detta exempel har fler lösningar än en, men för att inte exemplet ska bli för
komplicerat
utelämnas
dessa.
Genom
att
variera
parametrarna
för
poängtransformationsmatrisen och kostnaden för luckor kan andra optimala
stigar uppstå. Det gäller att följa vissa normer vid variation av parametrarna
för att en så representativ sekvensinpassning som möjligt ska kunna erhållas.
Ofta räcker det inte att ge samma kostnad för varje lucka. Enligt Korf m.fl.
[KYB03 s. 51] har det visat sig att raderingar och insättningar av nukleotider
tämligen ofta sker flera i rad . Detta betyder att om en lucka hittas i en viss
position i strängen, är sannolikheten stor att också den där på följande
position motsvarar en lucka, och så vidare. Needleman-Wunsch algoritmen
12
använder endast uniforma kostnader för luckor, vilket betyder att sekvenser
av luckor i en sträng inte prioriteras över enskilt utspridda luckor.
För att lösa detta problem har en metod som straffar långa sekvenser av
luckor lindrigare tillämpats med hjälp av en konkav straffpoängskala
[Mnei02a s. 4 - 5] samt [Mnei02b s. 1]. Metoden går i stort sett ut på att för
varje lucka i en sekvens av luckor, använda en straffkostnadsfunktion, γ (n) ,
som givet parametern av längden av lucksekvensen n bör satisfiera formel
3.1:
γ (n + 1) − γ (n) ≤ γ (n) − γ (n − 1)
Formel 3.1. Se t.ex [Mnei02b s. 1].
Denna metod har dock en körtid O( MN 2 + M 2 N ) (där M och N motsvarar
nukleotidsekvensernas längder) eftersom för varje position i matrisen bör en
extra kolumn och rad itereras för kontroll av luckor, vilket medför att
metoden baserad på konkav straffpoängskala är mycket långsammare än den
enligt Needleman-Wunsch presenterade algoritmen.
γ (n) =e + (n − 1)d där n ≥ 1
Formel 3.2. Se [Mnei02b s. 1].
För att snabba upp körtiden kan en affin straffkostnadsalgoritm användas.
Observera att kostnaden som returneras från γ (n) dras av från den totala
poängen som ska maximeras. Enligt denna metod straffas endast den första
luckan hårt och återstoden av de på varandra följande luckorna straffas
uniformt samt även lindrigare än den inledande luckan. Formel 3.2 beskriver
denna funktion, där e är öppningskostnaden och d är utvidgningskostnaden för
en sekvens av luckor av längd n. Jämfört med den konkava metoden, som för
varje position i matrisen måste iterera över en kolumn och en rad, räcker det
att endast granska den senaste cellen för varje kolumn respektive rad. På detta
sätt kan körtiden fås tillbaka ner till O(MN).
13
Figur 3.2. Bild tagen från [Mnei02b s. 1] och modifierats genom att insättning
av linje A.
I figur 3.2 jämförs straffkostnaderna för de olika straffunktionerna. Den
horisontala axeln representerar antalet konsekutiva luckor och den vertikala
axeln representerar straffkostnaden. För bestraffningen av luckor enligt
Needleman-Wunsch
algoritmen
erhålls
en
linje
A.
En
konkav
straffkostnadfsunktion ger en böjd kurva B medan den affina straffunktionens
kostnad representeras som en bruten linje C. Som ur bilden kan observeras,
representerar
den
affina
straffkostnadsfunktionen
verkligheten
(som
motsvarar den konkava kurvan) bättre än vad Needleman-Wunsch gör.
Närmare information om hur den affina straffkostnadsalgoritmen fungerar
beskrivs i kapitel 4.1.2. Notera dock att algoritmerna som används i POY
strävar efter att minimera kostnaden till skillnad från de hittills presenterade
algoritmerna som försöker maximera poängen för goda inpassningar.
14
4. POY – ett praktiskt hjälpmedel inom kladistiken
POYär ett kladistiskt verktyg och används i syftet att skapa kladogram som
återspeglar evolutionens gång på ett så objektivt sätt som möjligt [Whee06 s.
34]. Både genetiska data och morfologiska data lämpar sig som underlag för
datoriserade beräkningar och även om POY kan utnyttja båda typerna av data,
kommer endast algoritmer för behandling av genetiska data att tas upp i denna
avhandling.
I POY har många funktioner inkluderats för konstruktion av optimala
kladogram och tonvikten har lagts på en teknik som kallas för Dynamisk
homologi (Dynamic Homology, se [Whee01]). Principen bakom dynamisk
homologi är att försöka avfärda så kallade ad-hoc teorier som kan uppstå vid
manuell sekvensinpassning och uppbyggnad av kladogram. Dynamisk
homologi består av två nästlade processer – kladogramsökning följt av
strängoptimering. Strängoptimeringsfunktionerna används vid sökning efter
bästa möjliga inpassningar av DNA-sekvenser för ett givet kladogram och
räknar ut dess specifika minsta kostnad. Den andra gruppen består av
heuristiska sökfunktioner vars uppgift är att välja ett sådant kladogram som
ger en optimal (minimerad) kostnad. Sökfunktionerna väljer ut olika
hypotetiska
topologier
som
efter
hand
kan
testas
med
hjälp
av
strängoptimeringsfunktionerna.
4.1. Optimering av strängar
Det finns flera olika strängoptimeringsalgoritmer implementerade i POY.
Några av dessa är Direct Optimization, se [Whee06 s.47], Iterative-pass
Optimization, se [Whee03a], Fixed-states Optimization, se [Whee99], samt
Search-based Optimization, se [Whee03b]. Den möjligtvis viktigaste och mest
använda algoritmen är direkt optimering och inom ramen av detta arbete
kommer endast denna algoritm att analyseras. En viktig sak att notera är att
strängar av DNA representeras i POY samt i kommande pseudoalgoritmer
15
som bitmängder. Exempelvis sekvensen (A/C)(G/C)TT , där (A/C)
denoterar ” A eller C ”, består av fyra baser som i IUPAC-format (se bilaga 4)
skrivs som MSTT där M motsvarar (A/C) och S motsvarar (G/C) . I POY
används inte IUPAC-formatet för att beskriva obestämda nukleotider eftersom
luckor i direkt optimering behandlas på samma sätt som de övriga baserna.
Det finns heller ingen överenskommelse för namn på ett format där en lucka
motsvarar ett femte tillstånd och pga. detta införs i denna avhandling ett nytt
namn på ett sådant format – IUPAC2.
4.1.1. Direkt optimering
Direkt optimering, se [Whee06 s. 47], är en metod som strävar mot att skapa
hypotetiska förfäder (Hypothetical Taxonomic Units, HTU) på ett sådant sätt
att det slutgiltiga kladogrammet får en minimal kostnad. Detta utförs med
hjälp av tvådimensionella stränginpassningsalgoritmer.
Direkt optimering är uppdelad i flera algoritmfaser vilka i tur och ordning
tillämpas på ett specifikt kladogram C som har valts av en lämplig algoritm
för kladogramsökning. Direkt optimering räknar ut kostnaden av C vilket
sökalgoritmen antingen förkastar eller behåller. När C av en sökalgoritm
skickats till direkt optimering för analys, görs först en nedpassage (down
pass) på det valda kladogrammet och slutligen en uppassage (up pass). En
nedpassage innebär en rekursiv postordningstraversering (genomlöpning av
det valda kladogrammet). Operationerna startar med utgångspunkt från löven,
vilket går ut på att kladogrammet successivt traverseras ned mot roten. Av
detta följer, med referens till figur 4.1.1.1, att HTU1 blir resultatet av den
första operationen (mellan OTU1 och OTU2) i en nedpassage. Den sista av en
serie operationer i en nedpassage avslutas vid roten (HTU4 i figur 4.1.1.1).
16
Figur 4.1.1.1. Rekursiv postordningstraversering av ett kladogram.
funktion Nedpassage( n )
krav: n är en nod av ett kladogram, initialt roten
krav: n.vänster.tecken är vänstra barnets sekvens av längd k
krav: n.höger.tecken är högra barnets sekvens av längd l
om n inte är löv {
Nedpassage(n.vänster)
Nedpassage(n.höger)
(TM, kostnad) <- Strängjämförelse(n.vänster.tecken,
n.höger.tecken)
(inpassning1, inpassning2) <- Bakåtspårning(TM, k, l)
n.vänster.tecken <- inpassning1
n.höger.tecken<- inpassning2
n.tecken <- Bäst_HTU(n.vänster.tecken, n.höger.tecken)
n.kostnad <- n.vänster.kostnad + n.höger.kostnad + kostnad
}
Algoritm 4.1.1.1. Nedpassage.
Varje nod som inte är ett löv (OTU) kommer först, enligt algoritm 4.1.1.1, att
traverseras via vänster förgrening tills en lövnod nås. Eftersom lövet inte har
några förgreningar hoppar man tillbaka till dess förfader n som nu traverseras
via den högra förgreningen. Om höger förgrening också visar sig vara ett löv
återvänder man till förfadern vilken motsvarar en HTU, för att utföra en
sekvensinpassning av dess avkommors DNA-strängar. En sådan inpassning
görs för att den minsta kostnaden för jämförelse av två sekvenser ska kunna
erhållas på ett objektivt sätt (notera att uttrycket kortaste distans kan
genomgående användas i jämliket med minsta kostnad). Inpassningen
förklaras närmare i algoritm 4.1.1.2. samt algoritm 4.1.2.1. I samma procedur
erhålls också en bakåtspårningsmatris TM som innehåller information om
vilka riktningar den optimala stigen har. Att inpassa två DNA-sekvenser
17
innebär en sökning efter bästa möjliga matchning av nukleotider genom att
vid behov införa luckor i någondera sekvensen. Ju bättre ett par justerade
sekvenser är inpassade, desto lägre kostnad (kortare distans) får jämförelsen. I
POY används en modifierad version av Needleman-Wunsch algoritmen som
tillämpas vid global sekvensinpassning. I originalversionen av den nämnda
algoritmen används ett poängmaximeringssystem till skillnad från algoritmen
implementerad i POY. Problemet med poängmaximering är att det kan vara
svårt att se när ett par av strängar matchar varandra exakt. Används ett
minimeringssystem närmar sig kostnaden i detta fall alltid noll för en god
inpassning, medan resultat av maximering varierar beroende på längden av de
sekvenser som jämförs.
funktion Strängjämförelse( S, T )
krav: S och T är sekvenser med längderna m och n
krav: M är en tom kostnadsmatris av storlek (m+1) x (n+1)
krav: TM är en tom bakåtspårningsmatris
M[0][0] <- 0
för varje i = 1 till m {
M[i][0] <- M[i – 1][0] + σ(S i , lucka)
}
för varje i = 1 till n {
M[0][i] <- M[0][i - 1] + σ(T i , lucka)
}
för varje i = 1 till m {
för varje j = 1 till n {
kostnad1 <- M[i – 1][j – 1] + σ(S i , T j )
kostnad2 <- M[i – 1][j] + σ(T i , lucka)
kostnad3 <- M[i][j – 1] + σ(lucka, S j )
M[i][j] <- min(kostnad1, kostnad2, kostnad3)
TM[i][j] <- 0
om kostnad1 = M[i][j] {
TM[i][j] = TM[i][j] ∪ BACKA_DIAGONALT
}
om kostnad2 = M[i][j] {
TM[i][j] = TM[i][j] ∪ BACKA_VÄNSTER
}
om kostnad3 = M[i][j] {
TM[i][j] = TM[i][j] ∪ BACKA_UPP
}
}
}
returnera (M[i][j], TM)
Algoritm 4.1.1.2. Strängjämförelse.
Algoritm 4.1.1.2 fungerar i princip på samma sätt som Needleman-Wunsch
algoritmen, där den eller de minsta kostnaderna av kostnad1 , kostnad2 och
18
väljs.
kostnad3
Kostnaderna
räknas
ut
med
hjälp
av
en
kostnadsomvandlingsfunktion σ, som till skillnad från Needleman-Wunsch
algoritmens
motsvarande
poängomvandlingsfunktion
minimerar
goda
jämförelser. Kostnadsomvandlingsfunktionen representeras av en matris som
ställer upp alla kombinationerna av kostnadsomvandlingar för nukleotidbaser
inklusive luckor i ett tvådimensionellt plan. Figur 2.2.1 är exempel på en
möjlig kostnadsomvandlingsmatris. I POY tillåts användaren själv definiera
egna kostnadsomvandlingsmatriser. Enligt Press m.fl. [PTVF07 s. 869], bör
en sund kostnadsomvandlingsmatris (distansmatris) TCM uppfylla följande
krav:
TCM(i, j) ≥ 0
(positiv)
TCM(i, i) = 0
(ingen egenkostnad)
TCM(i, j) = TCM(j, i)
(symmetrisk)
TCM(i, k) ≤ TCM(i, j) + TCM (j, k)
(triangelolikheten bör gälla)
Figur 4.1.1.2. Exempel på en sund kostnadsomvandlingsmatris.
Det som är viktigt att notera är att sekvenser (strängar) använder IUPAC2formatet (se kapitel 4.1) vilket betyder att varje tecken i en sekvens kan
representera flera nukleotider samtidigt. På så vis kan varje tecken enligt
IUPAC2-formatet ses som en mängd med maximalt fem element (A/C/G/T/-),
varför mängdoperationer kan tillämpas i pseudoalgoritmerna. Likt detta, tar
kostnadsomvandlingsfunktionen σ
två tecken av IUPAC2-format som
argument. För att kunna avgöra kostnaden mellan dessa två tecken jämförs
varje element i respektive tecken och den billigaste kombinationen väljs.
T.ex.,
om
IUPAC2-tecknen
som
ges
som
argument
till
kostnadsomvandlingsfunktionen representerar (A/C/G) samt (C/G/-) ,
19
erhålls 9 jämförelser. Enligt matrisen i figur 2.2.1 skulle den billigaste
kombinationen, G mot G, producera kostnaden 0 som då blir resultatet av
kostnadsomvandlingsfunktionens jämförelse.
I algoritm 4.1.1.2 håller Pekarmatrisen TM reda på alla kombinationer av
optimala stigar som kan uppstå. När första delen av inpassningen gjorts
returneras kostnaden samt pekarmatrisen och vidare kommer funktionen för
nedpassage att kalla på en bakåtspårningsfunktion (algoritm 4.1.1.3), som
används för att spåra up den eller de bästa stigarna som uppstår vid
strängjämförelsen (algoritm 4.1.1.2).
funktion Bakåtspårning( TM, S, T )
krav: S och T är två strängar med längderna m och n
krav: TM är en bakåtspårningsmatris som sparats vid inpassning
av S och T
krav: inpassningS samt inpassningT är initialt två tomma
sekvenser
medan m > 0 eller n > 0 {
om BACKA_DIAGONALT ∈ TM[m][n]{
inpassningS <- S m + inpassningS
inpassningT <- T n + inpassningT
m <- m – 1
n <- n – 1
}
annars om BACKA_VÄNSTER ∈ TM[m][n] {
inpassningT <- lucka + inpassningT
m <- m – 1
}
annars {
(* backa uppåt *)
inpassningS <- lucka + inpassningS
n <- n – 1
}
}
returnera (inpassningS, inpassningT)
Algoritm 4.1.1.3. Bakåtspårning.
Bakåtspårningen startar i den sista cellen i matrisen, M[m][n] . Om startcellen
har en pekare till den diagonala cellen M[m-1][n-1] fogas baserna i S m samt
T n till inpassningS samt inpassningT och sedan sker ett hopp till den
diagonala cellen. Om inte M[m][n] har en diagonalpekare, undersöks det om
det finns en pekare till vänstra cellen (och sist till cellen ovanför). Om en
pekare finns, fogas en lucka på samma sätt som baserna fogats. Luckorna
justerar strängarna enligt pekarmatrisens riktlinjer. Även om det finns flera
20
optimala stigar kommer bakåtspårningsalgoritmen endast att välja en av
dessa. Märk väl att en sammanfogning av baserna har högsta prioritet.
När bakåtspårningen är klar returneras de två inpassade sekvenserna. Den
sista fasen i nedpassagen går ut på att försöka bestämma den slutgiltiga
sekvensen för noden n , som är en hypotetisk anfader (HTU) till n.vänster och
n.höger. HTU:ns sekvens väljs genom kombination av dess avkommors
inpassade sekvenser enligt algoritm 4.1.1.4. Målet är att bilda en så
verklighetstrogen HTU av de två avkommorna som möjligt för att inte
information ska försvinna vid fortsatt analys vid tillbakagång mot trädets rot.
Som ett exempel, kan HTU1 i figur 4.1.1.1 beteckna n och OTU1 samt OTU2
beteckna n.vänster respektive n.höger . Men även HTU4 har HTU1 och
HTU3 som sina avkommor och för dessa interna noder tillämpas precis
samma algoritmer.
funktion HTU_estimering ( S, T )
krav: S samt T är två förinpassade sekvenser med längderna m
krav: TMP är ett tecken av IUPAC2-format
krav: HTU blir den slutliga hypotetiska anfadern till S och T
för varje i = 1 till m {
TMP = S i ∩ T i
om TMP = φ eller TMP i = lucka {
HTU i = S i ∪ T i
}
annars() {
HTU i = TMP
}
}
Algoritm 4.1.1.3. Estimering av HTU.
Algoritm 4.1.1.3 beskriver hur den hypotetiska förfaderns DNA-sekvens
bestäms. Till algoritmen skickas två redan inpassade sekvenser som argument.
Sekvenserna är lika långa och itereras igenom position för position. Om
tecknen som jämförs inte har något gemensamt element, eller om endast en
lucka är gemensam, sammanslås tecknen. Exempelvis unionen av tecknen
(A/G/-) och (C/-) resulterar i (A/C/G/-) som blir det nya tecknet i
position
i
i sekvensen HTU . Om tecknen däremot har en eller flera
gemensamma nukleotider, kommer snittet av de två mängderna att sparas
21
tecknet i position
i
i sekvensen HTU. Till exempel snittet av tecknen
(A/G/T/-) och (A/G) resulterar i (A/G).
S: AA C
G
G C T CCCC - G C AGG
|| |
|
| | | |||| |
| | | |||
T: AA T (G/-)(G/-)C G CCCC T
T G T AGG
|| |
|
| | | |||| |
| | | |||
HTU: AA(C/T) G
G C(G/T)CCCC(T/-)(T/-)G(C/T)AGG
Exempel 4.1.1.1.
Exempel 4.1.1.1 illustrerar hur en estimerad förfaders sekvens, HTU, kan ser
ut efter en sammanslagning av sekvenserna S och T.
Figur 4.1.1.3. Ett slutgiltigt träd efter en nedpassage.
Efter att nedpassagen blivit färdig kan exempelvis ett träd som figur 4.1.1.3
illustrerar uppstå. Trädet har då traverserats i enlighet med figur 4.1.1.1.
En uppassage kommer till slut att utföras för att rätta till möjliga
inkonsekventa sekvenser som uppstått till följd av tidigare resultat från
nedpassagen. En uppassage är en rekursiv preordningstraversering, vilket
innebär att operationerna, i motsats till nedpassage, startar i roten och går upp
mot löven. Istället för att operera mellan två avkommor av en nod, som sker i
nedpassagen, sker operationerna mellan avkomma och förfader. Dessa
operationer bestämmer, på samma sätt som i nedpassage, de bästa kostnaderna
och passar in sekvenserna mellan avkomman och förfadern. Denna avhandling
går inte djupare in på uppassage, eftersom vetskapen om hur nedpassage
22
fungerar kan anses som tillräcklig för att ge en bild av vad som sker i direkt
optimering och för att kunna förstå huvudprincipen bakom POY.
4.1.2. Direkt optimering – inpassning med affin straffkostnad
I kapitel 4.1.1 redogörs för hur en nedpassage som använder NeedlemanWunsch algoritmen vid sekvensinpassning fungerar. Men som i kapitel 3
påpekas, beskrivs evolutionen bättre om affina straffkostnader används vid
sekvensinpassningen. I detta kapitel presenteras två pseudoalgoritmer som
ersätter
algoritm
4.1.1.2
samt
algoritm
4.1.1.3
i
den
ursprungliga
nedpassagen. De övriga algoritmerna ändras inte, förutom att algoritmen för
nedpassage
tar
strängjämförelsen
emot
tre
stycken
bakåtspårningsmatriser
istället
för
en
skickar
och
dessa
vidare
från
till
bakåtspårningsalgoritmen. I POY är dessa tre integrerade i en enda matris för
att minska på minnesanvändningen, men det är lättare att beskriva
pseudoalgoritmerna
om
tre
åtskiljda
bakåtspårningsmatriser
används.
Observera dock att algoritm 4.1.2.1 är en grov förenkling av den algoritm som
är implementerad i POY.
funktion Strängjämförelse( S, T )
krav: S och T är sekvenser av längderna m och n
krav: MD, MH och MV är tomma kostnadsmatriser av storlek (m+1)
x (n+1)
krav: TMD, TMH och TMV är tomma bakåtspårningsmatriser
krav: e och d motsvarar öppningskostnad samt utvidgningskostnad
för luckor
MD[0][0], MH[0][0], MD[0][0] <- 0
för varje i = 1 till m {
MD[i][0], MV[i][0] <- ∞
MH[i][0] <- e + (i – 1) * d
}
för varje i = 1 till n {
MD[0][i], MH[0][i] <- ∞
MV[0][i] <- e + (i – 1) * d
}
för varje i = 1 till m {
för varje j = 1 till n {
tmp1 <- MD[i-1][j-1] + σ(S i , T j )
tmp2 <- MH[i-1][j-1] + σ(S i , T j )
tmp3 <- MV[i-1][j-1] + σ(S i , T j )
MD[i][j] <- min(tmp1, tmp2, tmp3)
om MD[i][j] = tmp1 {TMD[i][j] <- TMD[i][j] ∪ TILL_TMD}
om MD[i][j] = tmp2 {TMD[i][j] <- TMD[i][j] ∪ TILL_TMH}
23
om MD[i][j] = tmp3 {TMD[i][j]
tmp1 <- MD[i-1][j] + d,
tmp2 <- MH[i-1][j] + e
MH[i][j] <- min(tmp1, tmp2)
om MH[i][j] = tmp1 {TMH[i][j]
om MH[i][j] = tmp2 {TMH[i][j]
tmp1 <- MD[i][j-1] + d
tmp2 <- MV[i][j-1] + e
MV[i][j] <- min(tmp1, tmp2 )
om MV[i][j] = tmp1 {TMV[i][j]
om MV[i][j] = tmp2 {TMV[i][j]
<- TMD[i][j]
∪ TILL_TMV}
<- TMH[i][j]
<- TMH[i][j]
∪ TILL_TMD}
∪ TILL_TMH}
<- TMV[i][j]
<- TMV[i][j]
∪ TILL_TMD}
∪ TILL_TMV}
}
}
returnera (min(MD[i][j], MV[i][j], MH[i][j]), TMD, TMH, TMV)
Algoritm 4.1.2.1. Affin Strängjämförelse.
I algoritm 4.1.2.1 används tre kostnadsmatriser istället för endast en. I dessa
matriser ( MD, MH och MV ) sparas endast specialiserade kostnader. Detta
innebär att t.ex. MD (diagonalkostnadsmatrisen) endast lagrar kostnader för
jämförelser av par av nukleotider ( S i och T j ), men inte för jämförelser av
nukleotider och luckor (kostnader för insättning av luckor). Genom att
undersöka samma elementposition (i, j) för de tre matriserna och för varje
matriselement lägga till en omvandlingskostnad σ(S i , T j ) , kan det (eller de)
alternativ som producerar lägst kostnad väljas och ett (eller flera)
pekarelement läggas till i position (i, j) i pekarmatrisen TMD . Om t.ex. MH[i1][j-1]
+
σ(S i , T j ) producerar lägst kostnad vid en jämförelse, läggs
pekarelementet TILL_TMH till pekarmängden i TMD[i][j], vilket informerar
pekaren som används vid bakåtspårningen att hoppa till TMH[i-1][j-1] vid
kommande iteration . På så sätt bildas den stig som representerar inpassningen
av de två sekvenserna.
I MH (horisontalkostnadsmatrisen) lagras endast kostnaderna för horisontala
luckor. Eftersom MH representerar kostnader för infogning av luckor i
horisontalled (dvs. i S ), skiljer sig straffkostnaden drastiskt för ett hopp från
TMD
till TMH jämfört med straffkostnaden för itererad infogning av luckor i
MH .
Ett hopp från TMD till TMH innebär att en lucka tillfogas sekvensen S , vars
föregående tillsatta element motsvarar en nukleotidbas. På detta sätt uppstår
en första lucköppning som bör straffas hårdare än fortsatt infogning av luckor
24
i MH . Ett hopp från TMD till TMH straffas med en öppningskostnad d , medan en
kontinuerlig infogning endast ”straffas” med en utvidgningskostnad e som bör
vara lägre än öppningskostnaden.
MV ( vertikalkostnadsmatrisen )
MH
och TMV fungerar i princip på samma sätt som
och TMH , med enda skillnaden att MV lagrar kostnaderna för vertikala
luckor.
När strängjämförelsen är klar returneras den bästa kostnaden för den bästa
kostnadsmatrisen samt alla tre pekarmatriser.
funktion Bakåtspårning( TMD, TMH, TMV, S, T )
krav: TMD, TMH och TMV är fyllda bakåtspårningsmatriser
krav: S och T är två strängar med längderna m och n
krav: inpassningS samt inpassningT är initialt två tomma
sekvenser
krav: p pekar initialt på den pekarmatris (TM*) som gav bästa
svaret i strängjämförelsen
medan m > 0 eller n > 0 {
om TILL_TMD ∈ p {
inpassningS <- S m + inpassningS
inpassningT <- T n + inpassningT
m <- m – 1
n <- n – 1
p <- TMD[m][n]
}
annars om TILL_TMH ∈ p {
inpassningT <- lucka + inpassningT
m <- m – 1
p <- TMH[m][n]
}
annars{ (* TILL_TMV ∈ p *)
inpassningS <- lucka + inpassningS
n <- n – 1
p <- TMV[m][n]
}
}
returnera (inpassningS, inpassningT)
Algoritm 4.1.2.2. Affin bakåtspårning.
Den affina bakåtspårningsalgoritmen (algoritm 4.1.2.2) påminner mycket om
algoritm 4.1.1.3, men eftersom den nu presenterade algoritmen använder tre
bakåtspårningsmatriser behövs en särskild pekare p som kan hoppa mellan
dem. Pekaren p startar i position (m, n) för den bakåtspårningsmatris vars
motsvarande kostnadsmatris’ element i nämnda position gav upphov till bästa
25
kostnad. De så kallade pekarelementen TILL_TMD , TILL_TMH och TILL_TMV
avgör till vilken matris pekaren p hoppar för varje iteration. Sekvenserna
inpassningS
och inpassningT inpassas enligt stigen pekaren väljer och på
så sätt bildas den slutgiltiga inpassningen som sedan returneras.
4.2. Kladogramsökning
Dynamisk homologi är en kombination av kladogramsökning samt sträng- och
teckenoptimering. Sökningsalgoritmen väljer ut ett kladogram och skickar det
sedan vidare till tecken- eller strängoptimeringsalgoritmerna som beräknar
dess kostnad. Om trädet är billigare att konstruera än det eller de hittills
billigaste kladogrammen sparas det, annars förkastas det.
Den säkraste metoden att hitta ett globalt optimalt kladogram är genom
testning av alla möjliga kombinationer av topologier som kan uppstå. För fyra
OTU:n erhålls två olika trädstrukturer med 15 unika topologier, se figur
4.1.2.1. I figuren motsvarar a , b , c samt d godtyckliga OTU:n som kan
kombineras på 15 olika sätt. Trädstruktur A kan ordnas på 12 olika sätt medan
trädstruktur B på endast 3 sätt. Observera att topologin för ett kladogram
enligt definition inte ändras vid rotering av en eller flera godtyckliga HTU:ns
förgreningar, eftersom sekvensinpassning är en kommutativ operation:
inpassning( sträng1, sträng2 ) ⇔ inpassning( sträng2, sträng1 ).
26
Figur 4.2.1. 15 topologier av 4 OTU:n.
Antalet unika, rotade kladogram kan räknas ut med hjälp av Formel 4.2.1 av
Felsenstein [Fels78], där n motsvarar antalet OTU:n. Problemet med denna
taktik är att sökrymden ökar exponentiellt i takt med antalet taxa som läggs
till. För 10 OTU:n bör 40 320 stycken kladogram analyseras och redan för 20
OTU:n finns det 6 402 373 705 728 000 unika kladogram. Dessutom bör man
minnas att varje kladogram ska testas för kostnad, vilket för n antal noder
(vars sekvenser inte är uppdelade) kräver n - 1 antal inpassningar för varje
unikt kladogram.
(2 n − 5)!
2 ( n − 2)!
n−2
Formel 4.2.1. Formel enligt Felsenstein [Fels78].
För att lösa problemet har olika heuristiska metoder implementerats. En
vanlig startmetod i POY är att först generera ett eller flera wagnerträd
(Wagner Trees se bl.a. [Farr70]) som utgångspunkt för att sedan fortsätta
analysera och manipulera det eller de resulterande kladogrammen genom att
applicera metoder som t.ex. Subtree pruning and regrafting [Whee06 s. 65 –
66]. Av alla heuristiska metoder behandlas endast Wagners algoritm i denna
avhandling, vilket bör vara tillräckligt för att ge en helhetsbild av POY.
27
4.2.1. Wagnerträd
Wagnerträd byggs med start från ett initialt kladogram bestående av två, ur en
mängd valda OTU:n. Genom att systematiskt infoga nya OTU:n för varje
möjlig plats i trädet uppstår nya träd vars kostnader examineras. Det eller de
träd som frambringar de minsta kostnaderna analyseras vidare på samma sätt
tills alla OTU:n blivit testade. De två initiala noderna väljs, beroende på
skript (script), oftast slumpmässigt ur mängden icke-examinerade OTU:n.
funktion Wagner( ϕ )
krav: ϕ är totala mängden OTU:n (T 1 ,...,T n )
krav: τ är ett initialt träd med två OTU:n
krav: W är en stack av klara wagnerträd, initialt tom
krav: TMP är en stack, innehåller initialt τ
krav: TMP2 är en stack, initialt tom
medan TMP ≠ φ {
C <- pop(TMP)
om C.otus = ϕ {
push(W, C)
}
annars {
bäst <- ∞
för varje OTU t ∉ C {
för varje båge (w, v) i C {
gör en kopia NC = (NV, NE) av C
avlägsna (w, v) ur NE
bifoga t samt en HTU-nod x till NE
bifoga (w, x), (x, v) samt (x, t) till NE
om NC.kostnad < bäst {
töm TMP2
push(TMP2, NC)
bäst <- NC.kostnad
}
annars om NC.cost = bäst { push(TMP2, NC) }
}
}
medan TMP2 ≠ φ {
push(TMP, pop(TMP2))
}
}
}
returnera W
Algoritm 4.2.1.1. Uppbyggnad av wagnerträd
TMP
samt TMP2 är temporära stackar som används som stöd då wagnerträd
byggs upp. Som utgångspunkt väljer algoritmen ett slumpmässigt kladogram
bestående av två lövnoder genom att ta det översta trädet C från stacken TMP .
28
Initialt finns endast trädet τ i TMP . En kopia ( NC ) av det valda trädet C görs
och varje nod som inte finns i NC testas för varje båge E i C (Observera att för
ett träd C = (E, V), representerar E mängden bågar för C och V mängden
noder). Figur 4.2.1.1 är en illustration över hur en ny OTU t , ur mängden ϕ C,
läggs in i ett träd. Observera att en ny HTU x också läggs till för att
upprätthålla kraven för ett fullständigt binärt träd (se definition, kapitel 2.2).
Figur 4.2.1.1. Insättning av en OTU.
Varje båge i trädet C examineras. Efter att t lagts till i NC uppstår ett nytt träd.
För det nya trädet testas kostnaden med hjälp av någon tecken- eller
strängoptimeringsmetod, t.ex. direkt optimering. Om kostnaden är bättre än
det hittills bästa trädet töms TMP2 och det nya trädet läggs på stacken. Om
trädet har samma kostnad som det eller de hittills bästa träden töms inte
stacken. På så vis kan flera wagnerträd uppstå.
Då alla bågar i C har examinerats av t väljs nästa OTU och samma procedur
upprepas. Då alla noder i ϕ -C testats förs slutligen det eller alla de bästa
träden sparade i TMP2 över på TMP och det är dags att återgå till början av
algoritmen och ta nästa träd från stacken.
Slutligen, då det inte finns flera OTU:n att addera till trädet, läggs trädet på
stacken W . Då TMP är tom avslutas algoritmen och stacken W , som innehåller de
bästa wagnerträden, returneras. På det returnerade trädet eller träden
29
appliceras ofta flera andra heuristiska algoritmer eftersom kostnaden av
wagnerträd ofta tenderar att resultera i lokala minima. Detta beror på att alltid
de bästa underträden väljs direkt, även om det kunde vara lönsamt att också
testa andra utgångspunkter.
4.3. Skript
Ett POY-skript består av en rad kommandon som styr vilka algoritmer och
data som ska användas. I POY 4 tillåts användaren definiera sina egna skript.
Inom ramen för detta arbete kommer endast några få skriptkommandon, samt
en inskränkt inblick av funktionaliteterna, att presenteras i tur och ordning.
Dessa kommandon har valts pga. att de återspeglar resultat av tester som
gjorts under optimeringsarbetets gång.
Kommandot read, se [VVBW s. 65 – 71], används för inläsning av datafiler. I
detta arbete används read endast vid inläsning av nukleotidsekvenser, även
om t.ex. morfologiska data eller fördefinierade träd kan läsas in. Behandling
av andra data än nukleotidsekvenser är utom räckvidd för detta arbete och
kommer inte att tas upp.
Kommandot build, se [VVBW s. 47 – 48], används för att bygga upp ett antal
kladogram med hjälp av wagners algoritm. Som parameter anges antalet träd
som ska genereras. För varje omgång kommer startpunkten att vara ett
slumpmässigt par av OTU:n. Allt för höga värden på build-parametern
tenderar att producera flera träd av samma topologi, om organismerna som
jämförs är få till antalet.
Kommandot transform, se [VVBW s. 97 – 106], anger vilka optimeringsalgoritmer som ska tillämpas och vilka kostnader samt kostnadsmatriser som
ska
användas.
gap_opening:2)
Till
exempel
transform((all,
tcm:”tcm.txt”),
instruerar alla sökalgoritmer att utnyttja direkt optimering
vid examinering av trädkostnad. Dessutom specificeras en av användaren
30
fördefinierad
gap_opening:2
poängomvandlingsmatris,
tcm.txt,
vid
strängoptimering.
anger att den affina inpassningsalgoritmen ska användas med
en öppningskostnad av 2. Om ingen öppningskostnad är definierad, används
den manipulerade versionen av Needleman-Wunsch istället.
Kommandot swap, se [VVBW s. 92 – 96], tillämpar olika heuristiska
algoritmer på ett givet träd som t.ex. genererats av wagners algoritm.
Algoritmer såsom TBR (Tree Bisection and Reconnection) och SPR(Subtree
Pruning and Regrafting) [Whee06 s. 67 – 71] utnyttjas för att överkomma
lokala optima som ofta uppstår vid konstruktion av kladogram med Wagners
algoritm. Ett underträd skärs av för en viss region och införs på annat håll i
trädet, för att sedan testa kostnaden av det nya trädet med en tecken- eller
strängoptimeringsalgoritm vald av kommandot transform.
Vid sökning lagras de bästa träden i minnet. Exempelvis build(200)
lagrar 200 träd i minnet, även om flera är identiska. Genom att använda
kommandot select, se [VVBW s. 47 – 48], väljs endast det eller de bästa
träden med unika topologier, beroende på argumenten givna.
Kommandot report, se [VVBW s. 73 – 81], används för att ge olika slags
rapporter av analyserna i olika format. Detta är en snabb process som inte
påverkar körtiden för ett skript.
Kommandot quit, se [VVBW s. 64], används för att avsluta programmet.
read ("data.fas")
transform ((all, tcm:"tcm.txt"))
build (5)
swap ()
select ()
report ("data", graphtrees)
quit()
Exempel 4.3.1. Ett testskript.
Exempel 4.3.1 motsvarar ett skript som läser in data från en fil data.fas,
förknippar omvandlingskostnaderna med omvandlingskostnadsmatrisen och
31
bygger upp fem Wagnerträd. För varje träd som byggs tillämpas också TBR
samt SPR pga. swap-kommandot. Efter att dessa fem träd har skapats väljs de
unika träd som har bäst kostnad och skrivs ut till en postscriptfil data.ps. Efter
det avslutas skriptet. För mer information om hur skripten fungerar, se
[VVBW].
32
5. Kodoptimering
Beräkningsintensiva program kräver mycket processorkraft och kan ta långa
tider att exekvera. Eftersom körtiden är dyrbar, speciellt vid företag som kör
beräkningstunga applikationer, är det av stor vikt att algoritmerna som nyttjas
i programmen är så effektiva som möjligt. Om en algoritm inte är tillräckligt
optimerad kan körtiden, jämfört med en optimerad version av algoritmen,
resultera i en enorm skillnad. Det gäller att ta ett programs totala körtid i
beaktande – om en representativ körning av programmet tar två sekunder, ter
sig en optimering onödig. Om användaren däremot måste vänta flera timmar,
dagar eller till och med månader, kan en optimering vara ändamålsenlig.
Utvecklare måste överväga lönsamhet av optimering och samtidigt ta
utvecklingskostnader, kundkrav och tid i beaktande för att finna det optimala
tillvägagångssättet.
En optimering på kodnivå innebär en omstrukturering av ett program utan att
algoritmkoncepten ändras. Slutresultat producerade av ett optimerat program
bör inte avvika från slutresultat producerade av programmets originalversion.
Optimering på kodnivå handlar oftast om hur koden kan omstruktureras så att
minnesanvändningen sköts så effektivt som möjligt och så att processorns
pipeline hålls fri från blockeringar samt tömningar till följd av onödiga
beroendekedjor eller felaktiga hoppförutsägelser.
Programutvecklare bör försäkra sig om att optimala algoritmer valts förrän en
optimering på kodnivå lönar sig att utföra. En algoritm med kvadratisk körtid,
med avseende på antal element, kan möjligtvis optimeras att bli hundratals
eller tusentals gånger snabbare om den kan ersättas med en algoritm med
logaritmisk körtid. I sådana fall påverkar sällan en kodoptimering av
programmet nämnvärt, procentuellt sett.
33
5.1. Moderna processorarkitekturer
För att kunna optimera ett program på kodnivå, är det bra att ha en generell
insikt i den moderna processorarkitekturen. Genom att känna till processorns
funktionaliteter kan man lättare avgöra vilka optimeringar som ger bra
resultat. I dag finns flera olika processorarkitekturer, men i denna avhandling
behandlas endast x86-arkitekturen, se [Casa07], vilket förhoppningsvis ska ge
en viss insikt om hur processorerna som baserar sig på denna arkitektur
fungerar.
Processorkapaciteten
har
ökat
enormt
de
senaste
decennierna
och
utvecklingen följer fortfarande Moores lag, se [Moor65]. I slutet av 80-talet
lanserades x486-processorn. Enligt en kort artikel från CPU-INFO.COM
[Cpui06] hade denna banbrytande processor nya innovationer som en enkel
fem-stegs pipeline och ett integrerat cacheminne (cache memory, se
[Smit82]). Från och med detta har pipelinetekniken och cacheminnen
förbättrats drastiskt.
5.1.1. Pipelining
Prestanda för en processor mäts ofta i antal bearbetade instruktioner per
tidsenhet, vilket är beroende av processorns klockfrekvens (angett i MHz eller
GHz). Tiden det tar att processera komplexa instruktioner kan ta många
klockcykler, varför processorkraft går förlorad pga. sysslolösa logiska enheter
i processorn. Detta problem kan lösas med hjälp av en teknik som kallas
pipelining.
Figur 5.1.1.1. Utan pipeline.
34
Pipelining innebär att instruktionerna indelas och opereras på i olika faser, så
att flera instruktioner kan processeras samtidigt men i olika fragment. Detta
innebär inte att bearbetningstiden för en specifik instruktion blir kortare,
tvärtom kan en enskild instruktion ta längre tid på sig innan den passerat
pipelinen.
Det
som
gör
pipelining
effektivt
är
det
faktum
att
genomströmningen av instruktionerna ökar.
Figur 5.1.1.2. Med pipeline.
Figur 5.1.1.1 beskriver ett enkelt scenario där en pipeline inte används.
Instruktionerna delas upp i delar, så att varje del kan processeras i skilda
logiska enheter för varje klockcykel. I detta mycket enkla fall finns fyra
logiska enheter: hämta (fetch), avkoda (decode), exekvera (execute) och skriv
(write). I hämta-enheten hämtas en ny instruktion som står på tur i
instruktionscacheminnet (instruction cache, se [Pas02 s. 5]). I avkoda-enheten
avkodas instruktionen för att få reda på operation samt möjliga operander.
Exempelvis för instruktionen add r1,r2 erhålls värden ur registren r1 och
r2, som skickas vidare till exekvera-enheten. I exekvera-delen, som består av
flera olika funktionella enheter, sker den önskade operationen mellan
operanderna. I fallet med add-operationen aktiveras den aritmetisk-logiska
enheten för heltalsoperationer och producerar resultat av operationen. I skrivenheten skrivs resultat tillbaka till något register eller möjligtvis direkt till
minnet. När instruktionen passerat dessa logiska enheter är det dags för nästa
instruktion, som processeras på liknande sätt.
Problemet med det nyss nämnda tillvägagångssättet är att slöseri med resurser
uppstår om endast en av de fyra enheterna tillåts operera åt gången. Genom
tillämpning av pipelining har man kunna kringgå detta problem och flera
instruktioner kan processeras samtidigt (i enlighet med figur 5.1.1.2). När
instruktion 1 passerar hämta-enheten och går in i dekoda-enheten stagnerar
35
inte pipelinen, utan instruktion 2 läses i sin tur in i hämta-enheten. De logiska
enheterna är synkroniserade enligt klockfrekvensen, vilket innebär att längden
på klockcyklerna inte får vara kortare än tiden det tar för den långsammaste
enheten att processera färdigt. Detta kan vara problematiskt om man vill höja
processorns klockfrekvens, varför man i nyare processorarkitekturer har
övergått till superpipelining [Pime s.4 – 5]. Superpipelining innebär att de
olika logiska enheterna delas in i flera underenheter för att göra stegen i
pipelinen kortare. Genom att göra stegen kortare och mer välbalanserade,
minskar
också
tidsåtgången
i
respektive
enhet,
vilket
innebär
att
klockfrekvensen kan höjas. Pentium 4 är ett exempel på en processor som
implementerat en djup pipeline. Olika versioner av Pentium 4 processorn
finns, men den version som har djupaste pipelinen är Pentium 4 Prescott med
31 steg. Enligt Glaskowsky har dock prestandan knappt förbättrats med en så
djup pipeline [Glas04].
5.1.2. Superskalär pipelining
För att ytterligare förbättra prestanda kan pipelinen göras superskalär, vilket
innebär att grupper av instruktioner kan exekveras parallellt.
36
Figur 5.1.2.1. Superskalär pipeline.
Eftersom somliga instruktioner är mer komplexa än andra, kräver de flera steg
i
pipelinen
varav
vissa
steg
dessutom
kräver
flera
klockcykler.
Heltalsoperationer processeras oftast snabbare än flyttalsoperationer, varför
de förstnämnda har färre antal steg och kan ”träda ur” pipelinen i ett tidigare
skede. Moderna processorer har dock möjlighet att använda olika SIMDextensioner
som
vektoriserar
flyttalsinstruktioner
så
att
flera
flyttalsoperationer kan ske parallellt [Inte06 s. 10 – 12]. De skript som
används i körningarna av POY medför inte nyttjande av flyttalsberäkningar,
varför närmare information om SIMD-extensioner och liknande enheter
utelämnas i denna avhandling.
Figur 5.1.2.1 är en schematisk förenkling av en superskalär pipeline. I hämtaenheten hämtas tre instruktioner vilka sedan vardera avkodas och skickas, om
möjligt, parallellt till skilda pipelines i exekvera-enheten. Idealiska är sådana
fall, där instruktionerna som följer varandra inte kräver samma funktionella
enhet. Men problem uppstår om flera instruktioner matas in utan först
granskas. Om instruktionerna är på väg till samma funktionella enhet minskar
parallellismen eftersom andra funktionella enheter står sysslolösa. Detta beror
förstås också på specifika implementationer av processorer, eftersom vissa
processorer kan ha flera funktionella enheter av samma typ och på så vis
parallellisera operationer av samma slag.
För superskalära pipelines har dock en schemaläggare för instruktioner (outof-order execution mechanism) implementerats för att öka parallellismen.
Denna skedulerare väljer ett antal (i detta fall tre) instruktioner, som står på
kö i instruktionsminnet, på ett sådant sätt att de alla passar i skilda
funktionella enheter. Detta resulterar i att somliga instruktioner kan bli
exekverade långt före andra, varför schemaläggaren bör se till att välja sådana
instruktioner som inte är beroende av varandra och inte heller är beroende av
instruktioner längre fram i kön. Detta kan underlättas med hjälp av
intelligenta kompilatorer, som kan ordna instruktionerna på förhand och
underlätta
arbetet
för
processorn.
Men
också
den
som
skriver
ett
37
datorprogram bör känna till detta faktum och på så vis försöka undvika långa
kedjor av beroenden mellan variabler.
Om grupper om tre instruktioner parallellt kan matas in i pipelinen har
processorn i bästa fall en genomströmning av tre instruktioner per klockcykel.
Detta är en markant förbättring jämfört med icke-skalära pipelines. Dagens
processorer har förenat superpipelining med superskalär pipelining och kallas
endast pipelining.
X86-arkitekturen baserar sig på CISC-instruktioner, men har i modernare
processorer övergått till en RISC-liknande arkitektur på grund av att CISCinstruktioner finns i mängder och massor vilket gör det svårt vid tillämpandet
av en pipeline [RPF03 s.65-66]. Genom att avkoda dessa makroinstruktioner
(CISC-instruktioner) till enklare (men fler) mikroinstruktioner har strukturen
förenklats och pipelinen blir lättare att implementera.
5.1.3. Instruktionsberoenden och Hoppförutsägelser
Ett problem till följd av pipelining kan uppstå då en instruktion är beroende
av en annan instruktion som tar lång tid att processera. Även om kompilator
och instruktionsschemaläggare lägger stor tonvikt på att försöka gruppera om
instruktionerna så att de inte är beroende av varandra, uppstår ändå situationer
då delar av pipelinen måste vänta på resultat.
x = y / z;
x2 = x + 7;
Exempel 5.1.3.1.
I exempel 5.1.3.1 är x2 beroende av resultatet för operationen y/z . Division
är en mycket långsam heltalsoperation och tar många klockcykler att utföra.
Eftersom den andra instruktionen måste vänta på att divisionen ska
processeras klart, införs så kallade bubblor (bubbles) i pipelinen, vilka
motsvarar instruktioner som inte utför någonting. Ju längre en pipeline är,
desto flera instruktioner strömmar genom den med påföljd att risken för flera
38
instruktioner som är beroende av varandra blir större, vilket kan medföra
bubblor som stagnerar pipelinen. I vissa fall kan beroendetiderna förkortas
genom utnyttjande av mellanliggande resultat. Om exempelvis en addition
utförts och resultatet redan finns tillgängligt som mellanhandsresultat men
ännu inte skrivits till ett register, kan ändå en instruktion som beror på detta
resultat använda det. På så sätt kan en del bubblor elimineras.
För fall liknande det i exempel 5.1.3.1 ser antagligen både kompilator och
schemaläggare beroendet mellan x och x2 och kan möjligtvis vänta med att
föra
in
x ’s
instruktion
i
pipelinen.
Men
om
också
de påföljande
instruktionerna är beroende av x eller x2 , så att en kedja av beroenden bildas,
kan det hända att komplexiteten blir så hög att det blir svårt att undgå
införandet av bubblor.
Ett annat stort problem med pipelining uppstår när koden som processeras
innehåller många hoppinstruktioner. Ett hopp innebär att instruktioner från ett
annat ställe i koden bör hämtas och matas in i pipelinen. Om hoppet är
ovillkorligt (unconditional jump) kan processorn läsa in koden i förväg. Men
ofta är hoppen villkorliga (conditional jump), vilket leder till att processorn
inte på förhand kan veta vilken del av koden som ska läsas in. Exempel
5.1.3.2 illustrerar detta problem:
if (x > y) a = b;
Kompileras till liknande kod:
jämför x och y
om x > y, hoppa till L1
annars a = b
cmp x, y
jl L1
mov b, a
L1:. . . .
Exempel 5.1.3.1.
I exemplet visas hur koden kan se ut för en villkorlig sats. Stegen från att en
hoppinstruktion läses in till att den avkodas och exekveras kan vara många i
en djup pipeline. På så sätt, när instruktionen jl L1 kommer in i pipelinen
kan det ta flera klockcykler innan resultatet kan avläsas. Om inga andra
instruktioner
matas
in
i
pipelinen
under
tiden
den
villkorliga
39
hoppinstruktionen processeras, förloras mycket resurser. Enligt Patterson
[Pat03] består ca. 15% av ett programs totala antal instruktioner av villkorliga
hopp, varför det är viktigt att inte stagnera pipelinen för varje ovillkorligt
hopp. Dessbättre har olika hoppförutsägelsemetoder (branch predicition
methods) tillämpats för att minska på problemets omfattning. Genom att låta
processorn förutsäga, med hjälp av viss logik, kan den med varierade resultat
förutspå om ett hopp kommer att ske eller inte. På så sätt kan de påföljande
instruktionerna matas in som förr och processorn väljer den kod som mest
sannolikt står på kö. Om processorn gissar fel, som ofta händer då hoppen
sker enligt arbiträra mönster, måste pipelinen tömmas och register möjligtvis
återställas (återställning av registren beror på hur schemaläggaren matat in
instruktioner som egentligen borde ske efter hoppinstruktionen).
Att tömma pipelinen och starta på nytt tar mycket resurser. Därför har många
olika hoppförutsägelsemetoder tillämpats. De enklaste metoderna är former av
statiska hoppförutsägelser (static branch prediction), som t.ex. kan vara en
överrenskommelse att villkorliga hopp aldrig kommer att ske (branch never
taken). För slingor liknande den i exempel 5.1.3.3 och som inte innehåller
andra villkorliga satser än den i början, gissar "aldrig-hopp-tekniken" alltid
rätt (förutom sista iterationen). Detta beror på att programräknaren (program
counter) endast ”faller igenom” om villkoret är sant och i slutet av slingan
utförs ett ovillkorligt hopp tillbaka till början.
while (i < j){
x[i] = y[j];
i++;
j--;
}
Exempel 5.1.3.3.
Men metoden är väldigt primitiv och i värsta fall gissar den alltid fel, vilket är
fallet om den villkorliga satsen avslutar slingan istället för att påbörja den,
som i exempel 5.1.3.4. Observera dock att dessa resultat kan variera,
beroende på hur kompilatorer implementerar hoppinstruktionerna.
do {
40
x[i] = y[j];
i++;
j--;
} while (i < j)
Exempel 5.1.3.4.
Statisk hoppförutsägelses kan också förbättras – de hopp som går bakåt i
koden utförs alltid medan de hopp som går framåt aldrig utförs. Statisk
hoppförutsägning fungerar ofta som reservutväg för mer
avancerade
dynamiska metoder som används i dagens processorer.
Genom att använda metoder som utnyttjar mönsterigenkänning av de senaste
hoppen kan
ytterligare felgissningar elimineras. En två-nivå adaptiv
hoppförutsägelseenhet (two-level adaptive branch predictor) som beskrivs av
Fog [Fog07b s. 8 – 10] är en dynamisk mekanism som utnyttjar mönster i
saturerande räknare (saturation counters) vilka baserar sig på de senast
utförda hoppen sparade i en BTB-buffert (Branch Target Buffer), se
[LeSm84].
Figur 5.1.3.2. Saturerande räknare.
En saturerande räknare består av två bitar och kan likt figur 5.1.3.2
representera fyra möjliga lägen för gissning om ett hopp ska utföras eller inte:
starkt tagen (strongly taken), svagt tagen (weakly taken), svagt inte tagen
(weakly not taken) och starkt inte tagen (strongly not taken). Om en
saturerande räknare motsvarar läge starkt tagen eller svagt tagen, är
förutsägelsen att ett hopp utförs. Däremot om räknaren motsvarar läge svagt
inte tagen eller starkt inte tagen, är förutsägelsen att inte utföra hoppet. En
saturerande räknare fungerar som en Turingmaskin, se [LeWi00], och kan
endast ”befinna” sig i ett läge åt gången.
41
Ett exempel: om den saturerande räknaren motsvarar läge svagt tagen,
kommer processorn att förutspå ett hopp. Eftersom hoppet endast är en
gissning, kan det visa sig att fel beslut fattas. I sådana fall bör pipelinen
tömmas och räknaren uppdateras enligt ett ej taget hopp (ET, genom att ta ett
steg till höger enligt bilden). Den saturerande räknarens nya läge blir då svagt
inte tagen. Om hoppet förutspåddes rätt, kommer räknarens nya läge att
uppdateras till starkt tagen.
Figur 5.1.3.3. Två-nivå adaptiv hoppförutsägelsemekanism.
En
två-nivå
innehållande
adaptiv
2n
hoppförutsägelsemekanism
består
av
en
tabell
saturerande räknare, illustrerad i figur 5.1.3.3. Rätt
saturerande räknare väljs genom avkodning av mönstret i hopphistoriken (där
ettorna representerar hopp som utförts och nollorna representerar icke utförda
hopp). På så sätt finns en individuell saturerande räknare för varje möjligt
mönster som kan bildas av de n senaste hoppen. Efter en ”uppvärmning” kan
en två-nivå adaptiv hoppförutsägelsemekanism förutspå ganska avancerade
mönster som upprepar sig.
Om exempelvis en två-nivå adaptiv hoppförutsägelsemekanism används med
n = 2 och hoppen som utförs följer mönstret 0011-0011-0011-...-0011, kan
hoppen efter en kort inlärning förutsägas perfekt. 00 uppdaterar räknare
nummer noll enligt kommande hopp som i mönstret är 1. Uppdateringen gör
att den saturerande räknaren byter läge enligt ett taget hopp (T). 01
uppdaterar på samma sätt räknare nummer ett att flytta sig enligt taget hopp.
42
11 uppdaterar räknare tre att uppdateras enligt inte taget hopp och samma
gäller 10 som uppdaterar räknare nummer två. Efter denna första inlärning
kommer nu de saturerande räknarna att kunna förutspå hoppen för det nämnda
mönstret.
Enligt Fog [Fog07b s. 9] använder sig Pentium MMX, Pentium pro, Pentium 2
och Pentium 3 alla av en två-nivå adaptiv hoppförutsägelsemekanism med n =
4. Vidare påpekar Fog [Fog07b s. 12 – 13] att i nyare processorarkitekturer
har olika komplicerade neurala hoppförutsägelsemetoder och hybridversioner
implementerats. Information om sådana kommer dock inte att tas upp i denna
avhandling.
5.1.4. Cacheminnen
Enligt Moores lag har antalet transistorer i en integrerad krets fördubblats
vartannat år under flera decennier. Problemet är att utvecklingen av
primärminnet
(RAM-memory)
har
varit
och
fortfarande
är
märkbart
långsammare. Tiden det tar att flytta data mellan register och primärminnet är
en flaskhals för processorn. Åtkomsttiden för primärminnet är avsevärt längre
än åtkomsttiden för registren. Detta beror på att primärminnet är dynamiskt,
dvs. uppbyggt av kondensatorer vilka tar tid att ladda och inte kan tillträdas
lika ofta som statiska minnen. Fortsättningsvis är primärminnet beläget långt
utanför
processorkretsen,
vilket
bidrar
till
ytterligare
latens
för
minnesåtkomst.
För att minska tiden det tar att hämta och skriva data till primärminnet har
olika typer av cacheminnen implementerats. Ett cacheminne är ett litet minne
av statisk typ som oftast är integrerat i processorkärnan, vilket gör det mycket
effektivt. Cacheminnet fungerar som en buffert mellan processor och
primärminne på så sätt att de primärminnesområden som tillträds av
processorn kopieras till cacheminnet, om de inte redan finns där. När
processorn gör ett nytt tillträde till ett minnesområde kommer först
43
cacheminnet att granskas. Om minnesområdet, som processorn är i färd med
att tillträda, redan finns i cacheminnet kan data returneras direkt från eller
skrivas direkt till cacheminnet, vilket är mycket snabbare än ett tillträde till
primärminnet.
Det som gör cacheminnen så effektiva är att tekniken utnyttjar det faktum att
majoriteten
av
program
uppvisar
både
temporala
samt
spatiala
lokalitetsegenskaper. Spatial lokalitet betyder att om data i en viss
minnesposition tillträds, kommer också närliggande minnespositioner högst
sannolikt snart att tillträdas. I praktiken innebär detta att program ofta
tillträder data sekventiellt, vilket cacheminnen drar nytta av. Istället för att
endast returnera data för den minnescell som processorn efterfrågade,
returnerar primärminnet ett helt block av sekventiella data som sedan lagras i
en cache-rad (cache line). Storleken på en cache-rad kan variera mellan
processorer, men är ofta runt 64 eller 128 byte. Detta görs eftersom de
närliggande datacellerna, enligt principen om spatial lokalitet, sannolikt
kommer att tillträdas. Temporal lokalitet innebär att en minnesposition som
tillträtts vid en viss tidpunkt högst sannolikt kommer att tillträdas igen inom
kort. Dessa två lokalitetsegenskaper gör att cacheminnet kan hållas tämligen
litet och ändå täcka primärminnet till hög grad. Enligt Patterson [Patt03]
finner processorn de sökta minnespositionerna i cacheminnet för ca. 90 % av
minnestillträdena.
Det finns flera generella tekniker som cacheminnen baserar sig på. Den äldsta
tekniken, som inte längre används i dagens processorer, kallas för direkt
kartläggning (direct mapping, se [Pas02 s.9 – 12]). Enligt denna teknik
associeras varje cache-rad med blockadresser som matchar, enligt figur
5.1.3.3. Om cacheminnet är 64KiB stort och en cache-rad är 64 byte, kommer
cacheminnet att ha 1024 rader. Om på samma sätt primärminnet är 512 MiB
stort, finns det 8 388 608 stycken block, varav varje cache-rad kommer att
handhålla 8 192 stycken block. Detta innebär att en cache-rad n endast kan
kartlägga block som uppfyller kravet n = blocknummer mod (antal cacherader). Stora problem med denna metod uppstår för program som tillträder
små delar av primärminnet med hopplängder motsvarande multiplar av
44
cacheminnets storlek. I sådana fall kommer det att ske en cachemiss för varje
minnesåtkomst, även om cacheminnet inte är fullt. En cachemiss innebär att
en minnesposition som söks upp inte finns buffrad i cacheminnet, varför ett
nytt minnesblock måste hämtas från primärminnet.
Figur 5.1.3.3. direkt kartläggning.
En bättre teknik, som ofta används i dagens processorer, kallas för n-vägs
associativ kartläggning (set associative mapping, se [Pas02 s.13 – 14]). Nvägs associativ kartläggning är en blandning av direkt kartläggning och fullt
associativ kartläggning (fully associative mapping, se [Pas02 s.12]) Om
cacheminnet använde sig av fullt associativ kartläggning, kunde vilket som
helst block ersätta vilken cache-rad som helst, men problemet är att
hårdvarulogiken blir för komplex att implementera. N-vägs associativ
kartläggning är därför en kompromiss av direkt kartläggning och fullt
associativ kartläggning. I denna teknik grupperas cache-raderna till små
mängder (cache line sets). För dessa mängder kan sedan primärminnesblocken
justeras associativt, vilket gör att primärminnesblocken inte längre är totalt
bundna till vissa cacherader. Enligt Torres är Athlon 64, Sempron samt
Opteron exempel på processorer som har implementerat 16-vägs associativ
cache, medan Pentium 4, Pentium D och Core 2 Duo är exempel på
processorer som har implementerat 8-vägs associativ cache [Torr07 s. 9].
45
Enligt
(art.
G.T.
hardw.sec.
s.3)
stödde
redan
386DX-arkitekturen
cacheminnen. Cacheminnet var då inte inbyggt i processorn, utan kunde
endast läggas till i moderkortet. I dagens läge har flera typer av cacheminnen
integrerats i processorbrickan. Det minsta och snabbaste cacheminnet kallas
för L1-cache. L1-cache ligger nära tillhands för de mest kritiska delar av
processorns pipeline och är inte mycket långsammare att tillträda än vad
registren är. L1-cachen har en typisk storlek mellan 32KiB och 256KiB, men
är oftast uppdelad i instruktionscache och datacache. Instruktionscache lagrar
instruktioner som ska exekveras, vilket gör att processorn snabbt kan plocka
upp de instruktioner som ligger i kö. Datacache lagrar som ordet beskriver
data (även om instruktioner egentligen också kan ses som data).
L2-cache är ett litet större cacheminne, som oftast också är integrerat i
processorbrickan. Denna typ av cache är långsammare än L1-typen och dess
storlek varierar, beroende på processor, vanligtvis mellan 256KiB och 4MiB.
L2-cachen fungerar som en buffert mellan L1-cache och primärminnet (eller i
vissa fall kan ytterligare en L3-cache finnas externt på moderkortet). Om en
L1-cachemiss sker, dvs. processorn inte hittar det den söker i L1-cachen,
kommer L2-cacheminnet att förfrågas som följande.
Orsaken till att cacheminnen inte görs större beror på att de är mycket dyrare
att tillverka samt att de tar mycket mer plats än dynamiska minnen.
Dynamiska minnen är uppbyggda av kondensatorer, vilket gör att de kan
konstrueras mycket mer sammanpressat än statiska cacheminnen som baserar
sig på vipptekniken [Torr07 s.2]. Eftersom utrymmet på processorkretsen är
begränsat, blir det problematiskt att tillföra större mängder cacheminnen.
5.2. Kodprofilering
I enlighet med Pareto-principen, se [FaHa s.3 – 7], spenderas största delen av
tiden vid körning av ett datorprogram endast i små delar av koden. Genom att
46
profilera programmen med profileringsverktyg kan man lokalisera de
kodpartier eller funktioner som tar mest tid att utföra. För många av
profileringsverktygen baserar sig mätningarna på en samplingsprocess. Figur
5.2.1 är utdata av en profilering som gjorts med profileringsverktyget gprof,
se [FeSt00] för mer ingående detaljer. För denna profilering har en
samplingsperiod på 0.01 sekunder använts, vilket innebär att för en funktion
f, med en total körtid på x sekunder, blir antalet samplingar N = x/0.01. Detta
betyder att funktioner med korta totalkörtider utsätts för ett mindre antal
samplingar, vilket gör att resultaten högst sannolikt blir mer felaktiga än för
funktioner med längre körtider. Om man utgår från att gprof utför N sampel
och ger som resultat olika funktioners relativa andelar av tiden: p1, p2,..., pn , så
n
att
∑p
n
= 1 gäller, där n motsvarar totala antalet funktioner och pn motsvarar
i =1
relativ andel av den totala tiden som samplats för en funktion, kan man med
hjälp av formeln σ p =
pi (1 − pi )
, räkna ut σ p som motsvarar standardfel för
N
proportion (standard error of a proportion, se [Bray57]). Detta medför att ju
fler samplingar (N) görs, desto mindre blir felestimeringen av körtiden för
respektive funktion pi . Om gprof använder en samplingsperiod på 0.01
sekunder och kör i 10 sekunder, utförs 1000 samplingar. Om då t.ex. en
funktion f tar upp 20 % av körtiden, blir felet σ p =
0.20(1 − 0.20)
≈ 0.01 , dvs.
1000
ca 1 %. Över en tio sekunders period tar funktionen upp två sekunder av
körtiden, vilket innebär ett fel på ca. 0,02 sekunder.
Figur 5.2.1. Profilering av script1.
47
Om en viss funktion utmärker sig extra mycket i resultatet av profileringen är
den en bra utgångspunkt för optimering. Figur 5.2.1 är ett utdrag av en
profilering efter att en icke optimerad version av POY kört script1, bilaga 3A.
Ungefär 72 % av körtiden spenderas i algn_fill_row, vilket gör den till den
överlägset mest använda funktionen. Om tiden som spenderats i denna
funktion teoretiskt kunde försummas, skulle programmets körtid snabbas upp
med en faktor på ungefär 3.5. Eftersom algn_fill_row tar upp en så stor del av
programmets körtid, ter det sig naturligt att undersöka möjligheter att
optimera just denna funktion. I kapitel 5.3 behandlas optimeringsmetoderna
som tillämpats.
Figur 5.2.2. Profilering av script2.
På samma sätt utmärker sig funktionen algn_fill_row_aff i en profilering av
POY körande script2, vilket kan ses i figur 5.2.2. algn_fill_row_aff tar upp
över 77 % av totala körtiden, vilket gör den till en klar kandidat för
optimering.
5.3. Kodoptimeringsmetoder
Efter att man försäkrat sig om att de bästa algoritmerna för kodoptimering
används kan man tillämpa olika kända optimeringsmetoder för att ytterligare
snabba upp exekvering av koden.
48
Denna avhandling behandlar endast optimering på kodnivå, vilket innebär att
de optimerade algoritmerna i POY inte ändrat karaktär. Kompilatorn sköter
kodoptimeringen till en stor del (om den är flaggad för optimering, vilket är
ett grundantagande) läs [FeSt00 s.65 – 94] för mer information. Beroende på
val av kompilator varierar hastigheten av det kompilerade programmet. Det
gäller att hålla koden så enkel som möjligt för att kompilatorn lättare ska
kunna göra en bra optimering. Ju bättre en kompilator är, desto svårare blir
det att åstadkomma bra resultat genom att optimera koden manuellt. Ett annat
problem vid manuell optimering är att man ofta inte känner till vilka
optimeringar kompilatorn kommer att applicera. Det kan ibland krävas en
undersökning
av
objektkoden
för
att
kunna
fastställa
vilka
optimeringsmetoder kompilatorn använder för att själv inte göra onödig
optimering. Objektkod kan exempelvis visualiseras med hjälp av verktyget
objdump [Gatl01 s. 11 – 14], vilket är inkluderat i de flesta distributioner av
Linux.
Inom ramen av detta arbete har endast Gnus C-kompilator gcc 4.1.2 [Stal03]
använts.
Alla
kompileringar
har
gjorts
på
den
teoretiskt
bästa
optimeringsnivån -O3 (av nivåerna -O0 inga optimeringar, -O1 optimera till
en del, -O2 optimera mer och -O3) för att resultaten av körningarna ska
korrelera. Att jämföra resultat av olika optimeringsnivåer samt metoder som
tillämpas av gcc är utanför räckvidden av denna avhandling.
5.3.1. Kodanalys
För mätningar av POY-körningar har skripten script1, script2 och script2b
använts (se bilaga 3). Alla de skript som används involverar direkt optimering
(se. kapitel 4.1.1). Script2 samt script2b betonar affina straffkostnader (se.
kapitel 3) medan script1 endast använder de enligt Needleman-Wunsch
definierade straffkostnaderna.
49
Enligt kapitel 5.2 utmärkte sig funktionerna algn_fill_row (se bilaga 1A) och
algn_fill_row_aff (se bilaga 2A) i profileringar av POY, vilket kräver närmare
analys av funktionerna.
Det som utförs i algn_fill_row (se bilaga 1A) är att, likt Algoritm 4.1.1.2, i
cell mm[i] sparas minsta kostnaden för någondera aa , bb eller cc , se figur
5.3.1.1. Bakåtspårningsmatrisen dm består av celler representerade av heltal
(ints). Dessa celler fungerar som pekarmängder för kommande bakåtspårning
(se algoritm 4.1.1.3). Eftersom en cell består av ett 32-bitars heltal, kan den i
teorin representera 32 distinkta bitpekare pga. att en bitpekare endast tar upp
en bit. För algn_fill_row krävs dock endast tre bitar för att kunna representera
alla stigar. ALIGN ( 0x1 ), INSERT ( 0x2 ) samt DELETE ( 0x4 ) är överrenskomna
konstanterna som representerar bitpekarna.
Figur 5.3.1.1. Kostnadsval.
För varje anrop till algn_fill_row itereras en ny rad för matriserna. I POY har
minneskraven
minskats
genom
att
endast
använda
två
räckor
för
representation av kostnadsmatrisen. Detta görs pga. att endast de två senaste
raderna används för temporära lagringar av resultat. Ändamålet med en
kostnadsmatris är inget annat än att för den sista cellen spara slutkostnaden av
en inpassning.
50
Figur 5.3.1.2. Svårigheter i den affina funktionen.
Funktionen algn_fill_row_aff bygger på algoritm 4.1.1.2 och är mer komplex
än algn_fill_row, även om principerna är tämligen lika. I algoritm 4.1.1.2 har
en hel del detaljer utelämnats för att hålla problemet på en någorlunda trivial
nivå.
För
direkt
optimering
med
affina
straffkostnader
används
tre
kostnadsmatriser i likhet med algoritm 4.1.1.2. Var och en av dessa matriser
representeras, på samma sätt som i algn_fill_row, av två rader. I algoritm
4.1.1.2 används tre distinkta pekarmatriser. I den verkliga implementeringen
av POY är dessa matriser sammanslagna till endast en pekarmatris, dm . Likt
algn_fill_row består varje cell i matrisen av en 32-bitars int som kan
representera 32 bitpekare. I den affina pekarmatrisen finns det fler pekare
pga. sammanslagningen: ALIGN ( 0x1 ), INSERT ( 0x2 ) samt DELETE ( 0x4 ).
Dessa är konstanter som representerar bitpekare för diagonalpekarmatrisen,
medan ALIGN_V ( 0x8 ) och DELETE_V ( 0x32 ) är konstanter som representerar
bitpekare för vertikalpekarmatrisen. ALIGN_H ( 0x64 ) och INSERT_H ( 0x128 ) är
motsvarande bitpekare för horisontalpekarmatrisen.
Svårigheterna som uppstår i den affina funktionen grundar sig på val av
stigar. Figur 5.3.1.2 är ett exempel på detta problem: En stigövergång från en
cell (med totalkostnaden 10) i position i+1 , j till en cell i position i+2,j+1 ,
får kostnaden 1 om transformeringskostnaden mellan nukleotidbaserna A och
C är 1. På så sätt har cell i+2,j+1 möjligtvis en minsta totalkostnad på 11.
Men också en luckutvidgning från position i+1,j+1 till i+2,j+1
ger
totalresultatet 11. Observera att luckutvidgning används här istället för en
lucköppning. En luckutvidgning betyder att en stig fortsätter följa samma
51
spår, antingen horisontellt eller vertikalt, genom fortsatt infogning av luckor.
Problemet som nu uppstår är att de två stigarna som leder till X resulterar i
olika
kostnader.
Stigen
( i+1,j;
i+2,j+1;
lucköppningskostnad för övergång från i+2,j+1
i+3,j+1 )
kräver
en
till i +3,j+1 eftersom
övergången mellan i+1,j; i+2,j+1 inte gav upphov till en horisontell lucka,
och på så vis måste en första öppningskostnad för horisontallucka infogas. En
lucköppningskostnad på 4 resulterar i en totalkostnad 15 för cell X . Den andra
stigen ( i+1,j+1;
i+2,j+1;
i+3,j+1 )
har redan infogat flera luckor
sekventiellt, vilket betyder att även nästa övergång, som sker mellan cellerna
i+2,j+1
och i+3,j+1 , endast har en utvidgningskostnad 1 vilket ger X en
totalkostnad på 12. Eftersom den förstnämnda stigen gav upphov till en
kostnad på 15 bör den inte markeras. Det gäller för algoritmen att reda ut
dessa problem så att den billigaste stigen verkligen väljs och att pekare till
mindre optimala stigar inte uppstår. På grund av liknande specialfall består
algn_fill_row_aff av många invecklade villkorssatser.
I både algn_fill_row_aff och algn_fill_row finns slingor. Eftersom koden
innanför en slinga itereras flera gånger än koden utanför är det av större vikt
att göra optimeringar på koden på insidan. Om exempelvis en slinga i
medeltal itereras 30 gånger, tillträds koden innanför slingan 30 gånger oftare
än den på utsidan. Till följd av detta blir en optimering innanför slingan 30
gånger mer märkbar än en liknande optimering på utsidan.
I kapitel 5.3.2 – 5.3.5 presenteras några optimeringsmetoder som har snabbat
upp exekveringen av POY. Det finns många optimeringsmetoder som kan
beprövas, men vissa är mer passande för ändamålet.
5.3.2. Optimerad cacheminnesanvändning
För att råka ut för så få cachemissar som möjligt gäller det att göra
konsekutiva minnestillträden. Om data som är lagrat i primärminnet inte
tillträds sekventiellt utan genom hopp mellan olika minnespositioner kommer
52
onödig data att laddas in till cacheminnet. Om man har en räcka av strukturer
som traverseras sekventiellt men endast tillträder något element för varje
struktur, bör man tänka om. Ett sådant fall kunde optimeras till att ha en räcka
av strukturer med pekare till räckor av element, så att elementen av samma
typ kan ordnas konsekutivt i minnet. Om antalet pekare för strukturen blir för
stort kunde möjligtvis räckan av strukturer ersättas med endast en struktur.
För att hålla cachemissarna nere är det viktigt att arrangera data som tillträds
ofta så nära varandra som möjligt.
I originalversionen av POY (build 1908) är datastrukturerna väl ordnade. De
för inpassning ämnade paren av sekvenser som lagras i strukturerna är
ordnade direkt efter varandra. Samma sak gäller raderna som representerar
kostnadsmatriser. För algn_fill_row_aff finns tre kostnadsmatriser och för
effektivitets skull är dessa matriser representerade av två rader var. Eftersom
minneshanteringen sköts tämligen bra för element som tillträds i de mest
intensiva delarna av koden, kvarlämnas inte mycket att optimera.
En smärre optimering av algn_fill_row gav dock en aning bättre resultat,
genom byte av raderna sex till åtta för koden i bilaga 1A.
cc = pm[i] + c;
bb = mm[i - 1] + gap_row[i];
aa = pm[i - 1] + alg_row[i];
ändras till
aa = pm[i - 1] + alg_row[i];
bb = mm[i - 1] + gap_row[i];
cc = pm[i] + c;
Optimering 5.3.2.1.
I detta fall har en iteration för räckan pm gjorts mera konsekutiv. Eftersom
iterationerna framskrider har data för pm, enligt den gamla versionen,
tillträtts i ordningen:
..., pm[n-1], pm[n–1], pm[n + 1], pm[n],...
jämfört med den ändrade versionen:
..., pm[n-1], pm[n], pm[n], pm[n + 1],...
53
För den första versionen kan extra cachemissar ske för de gånger pm[n]
överskrider en cache-rad, men direkt efter överskrivs av den gamla raden pga.
att pm[n-1] laddas in och ersätter den gamla, för att igen efter ersättas av
pm[n].
5.3.3. Eliminering av villkorliga hopp
Som i kapitel 5.1.3 togs upp kan förutsägelser av villkorliga hopp gå fel,
vilket leder till att pipelinen måste tömmas. För en stram slinga med
villkorliga satser som resulterar i enkla hoppmönster kan processorn ofta
förutsäga hoppen. Om hoppen sker slumpmässigt, kommer processorn att
förutsäga fel cirka 50 % av fallen.
För processorarkitekturer från och med P6 har en ny instruktionstyp
introducerats, CMOV (conditional moves, se [Inte06 kap. 7.3.1.1]). CMOVinstruktioner fungerar på liknande sätt som en hoppinstruktion, med den
skillnaden att istället för att utföra ett hopp flyttas data i ett register till ett
annat (endast om villkoret för instruktionen är sant). Genom att ersätta
villkorliga hoppinstruktioner med CMOV:s kan hopp elimineras, vilket
exempel 5.3.3.1 illustrerar.
if (x > y) a = b;
Kompileras till liknande kod utan CMOV:s:
cmp x, y
jl L1
mov b, a
L1:. . . .
jämför x och y
om y < x, hoppa till L1
annars a = b
Med CMOV:s kan hoppet elimineras och antalet instruktioner minskas:
cmp x, y
movl b, a
jämför x och y
om y < x, lägg a = b
Exempel 5.3.3.1.
Ett problem med villkorsinstruktioner som CMOV och SETCC är att de
avläser statusflaggor som finns i ett EFLAGS-register [Inte06 kap. 3.4.3]
54
(SETCC är en instruktionstyp som fungerar på samma sätt som CMOVinstruktioner, med skillnaden att destinationsregistret läggs till värdet 0x1 om
villkoret för instruktionen är sann). Resultat av jämförelseoperationer (t.ex.
cmp) lagras i dessa register, men även resultat baserade på andra
instruktioner. Det gäller att hålla noga reda på vilka instruktioner som rör
dessa flaggor för att inte felaktigheter ska ske. I exempel 5.3.3.1 ändrar
additionsinstruktionen på EFLAGS-registret, vilket gör att cmovle inte kan
använda utgångsresultatet av jämförelsen på första raden. Till följd av detta
blir det svårt att införa allt för många instruktioner mellan en ämnad
instruktion och en CMOV-instruktion.
cmp r1, r2
add r3, r4
cmovle r4, r5
Exempel 5.3.3.1.
Ett annat problem med SETCC- och CMOV-instruktionerna är att de är endast
effektiva för enklare former av villkorssatser. Med hjälp av villkorliga
hoppinstruktioner kan hopp ske över långa kodblock, men att skapa samma
typ av kod som sekventiellt kan itereras utan villkorliga hopp kan rent av vara
omöjligt.
En idé kunde vara att införa nya typer villkorliga instruktioner, t.ex. ADDCC
som likt CMOV och SETCC skulle utföra en operation. Resultat av ADDCC
skulle inte ändra på statusflaggorna (eller möjligtvis endast ändra på vissa
statusflaggor), så att beroenden av statusflaggorna kunde minska.
För optimering av algn_fill_row samt algn_fill_row_aff har dock CMOV- och
SETCC-instruktioner implementeras med goda resultat. Det som utförs i
algn_fill_row (raderna 10-50 i bilaga 1A) är att det minsta värdet av aa , bb
och cc sparas i mm[i]. I dm[i] flaggas de bitar som motsvarar variablerna
innehållande det minsta värdet. Om t.ex. aa är 7, bb är 7 och cc är 9, sparas 7
i mm[i] och bitpositionerna som motsvarar aa och bb flaggas i dm[i ] , i
enlighet med figur 5.3.3.1. En flaggning innebär att bit i en position läggs till
55
ett. ALIGN och INSERT är flaggkonstanter för aa och bb. Eftersom både aa och
bb
innehöll lägsta värdet, får dm[i] värdet 0x3 (resultatet av 0x1 | 0x2 ) efter
en logisk OR-operation dm[i] = (ALIGN | INSERT);.
Figur 5.3.3.1. Bitarna ordnade i ett heltalselement dm[i].
algn_fill_row innehåller många if-satser, vilka efter kompilering med gcc
(version 4.1.2) endast producerar villkorliga hopp. Genom att använda inline
assembly, se [Cole], kan man själv införa instruktioner på låg nivå, vilket har
gjorts i optimering 5.3.3.1:
__asm__(
"cmp %0, %1\n\t"
"cmovg %0, %1\n\t"
"mov $0x0, %0\n\t"
"cmovle %4, %0\n\t"
"setge %b0\n\t"
"cmp %1, %2\n\t"
"cmovl %2, %1\n\t"
"mov $0x4, %2\n\t"
"cmovl %3, %0\n\t"
"cmovg %3, %2\n\t"
"add %b0, %h0\n\t"
"add %h0, %b2\n\t"
: "+Q" (aa), "+r" (bb), "+Q" (cc)
: "r" (ZERO), "r" (TWO)
);
//
//
//
//
//
//
//
//
//
//
//
//
//
//
mm[i] = bb;
dm[i] = cc;
// rad 15
// rad 16
rad
rad
rad
rad
rad
rad
rad
rad
rad
rad
rad
rad
rad
rad
1
2
3
4
5
6
7
8
9
10
11
12
13
14
Optimering 5.3.3.1.
Optimering 5.3.3.1 har i princip ersatt raderna 10-50 i bilaga 1A. Det som
också bör tilläggas är att dessa konstanter också tillförts början av funktionen:
register const int TWO = 0x200;
register const int ZERO = 0;
I optimering 5.3.3.1 har alla villkorliga hopp eliminerats till följd av ifsatserna och programmet fungerar enligt följande:
56
Rad 13:
Följt av första kolontecknet ”:” kommer utdataregistren. Dessa register
fungerar både som käll- samt destinationsregister. Variabeln aa binds till det
dynamiska registret %0 , bb till register %1 och cc binds till register %2 . Exakta
register tilldelas inte för att ge kompilatorn frihet att välja de mest lämpliga
registren för variablerna då objektkoden genereras.
Rad 14:
Följt av nästa kolontecken kommer indataregistren. I dessa register binds
konstanterna TWO med register %3 och ZERO med register %4 .
Rad 1 och 2:
Jämför register %0 och %1 , vilka är associerade med variablerna aa och bb.
Om bb > aa , kopiera aa till register %1 . Detta betyder att bb skrivs över om
aa
är mindre. Vi vet nu att bb innehåller det minsta värdet av mängden ( aa ,
bb ).
Rad 3:
Register %0 (variabel aa) är nu onödigt och töms (ställs till noll).
Rad 4:
Om bb <= aa gäller, baserat på jämförelsen på rad 1, kopieras värdet för
konstanten TWO till aa . Värdet på konstanten TWO är 512, vilket är 0x200
hexadecimalt. I praktiken flaggas den tionde biten i registret sett från höger
(från den minst signifikanta delen av registret). Att den tionde biten flaggas
bör ses som att det högre av de två 8-bitars registren får värdet 2, vilket
motsvarar konstanten INSERT . Figur 5.3.3.1 beskriver delregistren av register
RAX (vilket valts godtyckligt för exemplet) för 64-bitars processorer (för 32bitars
processorarkitekturer
finns
inte
64-bitars
register,
RAX).
Om
konstanten TWO (0x200) läggs i register RAX, kommer register AL att vara 0
medan AH blir 2, vilket beror på att registren är sammansatta. Det går inte att
använda register EAX på samma gång som AX är i bruk, eftersom de
överlappar varandra. AH och AL är dock det enda par av register som inte
57
överlappar varandra, och på så vis kan man läsa från och skriva till dem utan
att information försvinner. Problemet är att inte alla register stöder
byteregistren. Inte heller alla instruktioner, som t.ex. CMOV, kan skriva till
det högre av de två byteregistren. Därför utförs skrivningen till AH genom att
bifoga åtta extra bitar till den minst signifikanta delen.
Figur 5.3.3.1. Register RAX (EAX).
Rad 5:
Om bb >= aa för jämförelsen på rad 1 gäller, flagga den minst signifikanta
biten i det låga byteregistret %b0 . På detta sätt kan vi utnyttja båda
byteregistren för register %0 . Om den minst signifikanta biten läggs till 1
betyder det att register %b0 får värdet motsvarande konstanten ALIGN .
Rad 6 och 7:
Jämför cc och bb . Om cc < bb , kopiera cc till bb :s register. Vi vet nu att bb
innehåller det minsta värdet, min( aa , bb , cc ), av de tre ursprungliga
variablerna.
Rad 8
Vi kan nu lägga värdet 4 (0x4 hexadecimalt) i register %2 , eftersom cc inte är
till mera nytta. Observera att 4 motsvarar värdet på konstanten DELETE , se
figur 5.3.3.1. Genom att tilldela cc detta värde är tredje biten i registret, sett
från höger, flaggad.
Rad 9
Om cc < bb gäller för jämförelsen på rad 6, töms registret (alla bitar läggs
till noll). Detta utförs genom att flytta register %3 , som motsvarar konstanten
58
ZERO ,
till register %0 (som är bundet till variabel aa ). Eftersom CMOV-
instruktioner endast kan kopiera register till register, måste konstanten ZERO
bindas till ett register.
Rad 10
Om cc < bb , nollställ register %2 ( cc ).
Rad 11
Addera register %0b och %0h och spara resultatet i det höga byteregistret %0b .
Rad 12
Addera det höga byteregistret %0h med det låga byteregistret %2b och spara
resultatet i det låga byteregistret %2b . Detta innebär att variabel cc innehåller
bitmönstret som beskriver bitpekarna i en cell för pekarmatrisen dm .
Rad 15 och 16
Här tilldelas mm[i] värdet för bb , som nu är det garanterat minsta värdet av de
ursprungliga variablerna aa , bb och cc . Bitmönstret (som representerar
pekarmängden) som sparats i cc tillskrivs också en cell i dm-matrisen. Bilaga
1B motsvarar nu koden efter att optimering 5.3.2.1 och optimering 5.3.3.1
utförts.
Även algn_fill_row_aff har optimerats på detta sätt. I bilaga 2A har raderna
38 – 80 ersatts med optimering 5.3.3.1 (resultatet kan observeras i Bilaga 2B),
med enda skillnaden att rad 15 i optimeringen ändrats till:
mm[i] |= bb;
5.3.4. Upprullning av slingor
Upprullning av slingor (loop unrolling) är en optimeringsteknik som
mångfaldigar blocket innanför en slinga. Detta görs för att kunna minska på
antalet hoppinstruktioner och jämförelseinstruktioner relativt till det övriga
59
antalet instruktioner, vilket medför en minskning av totala antalet bearbetade
instruktioner.
Vidare
kan
upprullning
av
slingor
öka
på
instruktionsparallellismen genom att låta flera instruktioner, som tidigare
exekverades i skilda iterationer, exekveras parallellt [AMD04 s.151].
Speciellt snäva slingor är viktiga att rulla upp. exempel 5.3.4.1 är typiskt
exempel på hur en upprullning kan se ut. Även om koden blir längre (och
fulare) blir den ofta snabbare, beroende på antalet iterationer för slingan. Ju
fler iterationer, desto mindre barlast läggs till. När slingan rullas upp som i
exempel 5.3.4.1, bör ett avslutande fall som tar hand om de sista iterationerna
tilläggas.
Före upprullning:
while (i < j){
i++;
x[i] = a[i] + i;
}
Efter upprullning:
while (i < j - 3){
i += 4;
x[i-3] = a[i-3] + i
x[i-2] = a[i-2] + i
x[i-1] = a[i-1] + i
x[i] = a[i] + i
}
while (i < j){ // behandla återstoden
i++;
x[i] = a[i] + i;
}
Exempel 5.3.4.1.
Resultatet av upprullningen i exempel 5.3.4.1 är att slingan itereras fyra
gånger mer sällan. Det betyder att fyra gånger färre jämförelser evalueras och
fyra gånger färre hopp utförs. Även parallellismen ökar för denna slinga
eftersom operationerna är totalt oberoende av varandra. Exekveringsenheten
kan nu lättare behandla dessa instruktioner parallellt, beroende på hur många
heltalsexekveringsenheter processorn har. För en slinga liten som denna kan
60
farten snabbas upp flerfaldigt. Om däremot antalet iterationer för slingan är
lågt, kan den ta längre tid att köra än utan upprullning. Om t.ex.
variabelvärdet för j sällan är större än 4 kommer den första delen av
upprullningen inte att köras vilket medför överflödiga villkorsinstruktioner.
I exempel 5.3.4.1 har slingan endast rullats upp med en faktor 4. Ju högre
upprullningsfaktor slingan har, desto mindre barlast i form av extra hopp och
evalueringar läggs till. Problem vid allt för höga upprullningar är att
instruktionscacheminnet blir mer belastat och att instruktionscachemissar
lättare uppstår som följd. Om slingan är belägen i en inline-funktion, kan det
vara anledning att hålla upprullningsfaktorn på en lägre nivå, eftersom alltför
stora kodpartier inte blir ”inlinade”.
Slingan i funktionen algn_fill_row har effektivt rullats upp med en faktor 8,
vilket kan observeras i bilaga 1C (som är en vidareoptimering av koden i
bilaga 1B). En annan smärre optimering som utförts i samband med
upprullningen är att operationen bb = mm[i – 1] flyttats ut ur slingan. Att
tilldela värdet av mm[n] till bb inne i slingan är onödigt, eftersom det
mellanliggande resultatet av bb redan motsvarar mm[n], vilket figur 5.3.4.1
illustrerar. Genom att flytta ut bb = mm[i-1] , sker denna tilldelning endast en
gång per funktionsanrop.
Figur 5.3.4.1. Utnyttjande av bb som mellanliggande resultat.
Funktionen
algn_fill_row_aff
har
däremot
inte
optimerats
med
upprullningstekniken, eftersom koden är för komplex. Blocket i slingan är
alldeles för stort för att en upprullning ska ha en positiv inverkan.
61
5.3.5. Övriga optimeringar
För funktionen algn_fill_row_aff har det varit svårt att hitta någon lämplig
optimeringsteknik. I detta kapitel presenteras endast två mindre optimeringar.
Detta är den mest optimerade versionen av algn_fill_row_aff (bilaga 2C) och
är en fortsatt optimering av koden i bilaga 2B.
Den första av optimeringarna som gjorts är att raderna 7-18 i bilaga 2B
flyttats utanför slingan (se bilaga 2C) pga. att variablerna c och cprev inte
ändras under iterationerna. Eftersom variablerna inte ändras, är det överflöd
att testa villkoren för varje iteration. Även om processorn kan förutsäga
hoppmönstren, kan en del barlast elimineras genom att ur slingan flytta ut
kodpartier som alltid genererar samma resultat. Ett problem är att pdnmm[i]
och pm[i] beror på variabeln i som ändras inne i slingan. Detta har lösts
genom införsel av variablerna cc_val , aa_val samt dd_val . Dessa variabler
tilldelas värden endast en gång och kan sedan utnyttjas som om de vore
konstanter, vilket sker på raderna 29, 30 och 42, se bilaga 2C.
Den andra optimeringen baserar sig till en del på testning. Det har visat sig att
om datasträngarna som representerar DNA är lika långa kommer villkoret på
rad 27 i bilaga 2B alltid att gälla (Om t.ex. kostnadsomvandlingsmatrisen i
bilaga 3F används, gäller påståendet för varje iteration av slingan). Om
villkoret på rad 27 alltid gäller är det underförstått att villkoren på raderna 19
och 23 inte gör det. I sådana fall evalueras de i onödan pga. att de är belägna
först i den sammanhängande sekvensen av villkorssatser. Genom att istället
flytta det sista blocket först kommer det alltid att evalueras först. För att ta
reda på vilka villkorsuttryck som bör gälla för villkorssatsen som flyttas först,
kan en sanningstabell som tabell 5.3.5.1 bildas. Man kan se från tabellen att
blocket för villkoret på rad 27 i bilaga 2B utförs om gap_row[i] == 0 &&
gap_row[i-1] == 0
eller om gap_row[i] > 0 && gap_row[i-1] > 0
gäller (notera att gap_row[n] inte får innehålla negativa värden). Vi får då ett
villkor som ser ut som detta:
if ((gap_row[i] == 0 && gap_row[i-1] == 0) || (gap_row[i] > 0
&& gap_row[i-1] > 0)){}
62
En ytterligare optimering kan göras genom att dela upp denna villkorssats så
att en villkorssekvens på totalt fyra block bildas. Detta kan göras eftersom för
lika långa DNA-sekvenser är uttrycket gap_row[i] == 0 && gap_row[i-1]
== 0
aldrig sant, varför motsvarande villkorssats kan flyttas sist i sekvensen.
Tabell 5.3.5.1.
gap_row[i]
0
0
>0
>0
gap_row[i-1]
0
>0
0
>0
sista blocket körs
ja
nej
nej
ja
Till sist kan man dra nytta av att satsen som flyttades först uppdelats i två
satser genom att uppdela den problematiska villkorssatsen på rad 31 ( bilaga
2B), så att den kan införas i sekvensen. Man kan dessutom utnyttja det faktum
att variablerna c och cprev inte ändras i slingan, varför satsen kan för den
första villkorssatsen i sekvensen ersättas med aa += aa_val . Detta knep
snabbar upp exekveringen ytterligare om första satsen i sekvensen oftast
gäller. Den slutgiltiga optimeringen motsvarar raderna 39-62 i bilaga 2C. Om
DNA-sekvenserna för organismerna som jämförs är lika långa kommer endast
det första blocket av sekvensen motsvarande raderna 39-62 i bilaga 2C att
utföras. I sådana fall kan den problematiska villkorssatsen, som före
optimeringen fanns på rad 31 i bilaga 2B, anses eliminerad. Med problematisk
villkorssats avses här en sats med många villkorsuttryck. En sådan kan orsaka
konflikter i BTB-bufferten, se [AMD s.16].
Även för DNA-sekvenser vars längder skiljer sig tämligen till ordentligt
mycket, som t.ex. sekvenserna i bilaga 3G, visar det sig att denna optimering
fungerar bra (optimeringsresultaten kan avläsas i kapitel 6.2).
63
6. Resultat
I detta kapitel presenteras resultat och utdata. För tidsmätningarna har POY
körts med tre olika skript, vilka finns i bilaga 3. Alla tre skripten utnyttjar
konceptet med direkt optimering, se kapitel 4.1.1 och 4.1.2.
Script1 använder inte affina straffkostnader för luckor, vilket medför att
algn_fill_row är den funktion som testas. Användande script1 bearbetas
DNA-sekvenserna hämtade ur filen data_script1.fas, bilaga 3E.
Script2 använder affina straffkostnader, vilket betyder att funktionen
algn_fill_row_aff är funktionen ämnad för testning. Användande script2
bearbetas
DNA-sekvenserna
hämtade
ur
data_script2.fas.
Alla
DNA-
sekvenserna är lika långa.
Script2b fungerar som script2, men DNA-sekvenserna i data_script2b.fas är
kortare och varierar i längd.
6.1. Utdata
I detta kapitel presenteras de utdata som uppkommit till följd av POYkörningar utnyttjande skripten i bilaga 3. Resultaten som erhållits är grafiska
representationer
av
kladogram.
Organismerna
som
presenteras
i
kladogrammen är ordnade enligt graden av likhet baserade på DNAsekvenserna skripten indikerat. Notera att utdata inte har påverkats av
optimeringarna gjorda i POY.
64
Figur 6.1.1. Utdata för script1.
Figur 6.1.1 är resultatet av en körning av script1.
Figur 6.1.2. Utdata för script2.
Figur 6.1.1 är resultatet av en körning av script2. Utdata för körningen
motsvarar tre unika kladogram, alla med samma kostnad. Observera att
kladogrammen skiljer sig tämligen mycket från varandra, vilket kan bero på
dåliga indata eller möjligtvis på att skriptet är för enkelt.
65
Figur 6.1.3. Utdata för script2b.
Figur 6.1.1 är resultatet av en körning av script2b. Märk att resultatet skiljer
sig kraftigt från alla kladogram producerade användande script2. Detta kan
bero på att DNA-sekvenserna som användes var för korta, eller möjligtvis att
skriptet är för enkelt.
6.2. Optimeringsresultat
För att få en uppfattning om hur mycket snabbare POY har blivit för varje
optimering som utförts, presenteras körtiderna i jämförbara tabeller. Tiderna
som presenteras motsvarar effektiv processortid. Optimeringsresultaten
presenteras i samma ordningsföljd som de togs upp i kapitel 5.3, eftersom de
baserar sig på varandra. Observera att varje optimering bygger vidare på de
föregående presenterade optimeringarna.
Alla tester har utförts i Fedora Core 4 vilket är ett Linuxbaserat
operativsystem. Som underlag för körningarna har en Core 2 Duo processor
fungerat, för vilken kärnorna har en maximal klockningshastighet på 2,40
GHz. Processorn har 4 MiB integrerat L2-cache, vilket de två kärnorna delar.
Tabell 6.2.1. Originalversionen av POY, körtider i sekunder.
körning 1
körning 2
script1
53,666
53,745
script2
40,974
41,111
script2b
25,185
25,299
66
körning 3
medeltal
53,698
53,703
41,113
41,066
25,213
25,232
I tabell 6.2.1 presenteras körtiderna för originalversionen av POY. Script1 tar
längst tid att processera, eftersom antalet organismer att ordna är flest, se
bilaga 3F. För script2 och script2b skiljer sig tiderna främst pga. att DNAsekvenserna som används i samband med script2b är mycket kortare än de för
script2.
Tabell 6.2.2. Optimering 5.3.2.1, körtider i sekunder.
körning 1
körning 2
körning 3
medeltal
script1
52,822
52,744
52,761
52,776
De små ändringar som gjordes i kapitel 5.3.2 gav, efter exekvering av script1,
upphov till resultaten presenterade i tabell 6.2.2. Optimeringen snabbade upp
medeltalskörningen med en faktor på ca. 1,02 jämfört med originalversionen
av POY. Denna optimering applicerades endast på funktionen algn_fill_row.
(Notera att en uppsnabbning med en faktor x bör tolkas som att programmet
exekveras x gånger snabbare.)
Tabell 6.2.3. Optimering 5.3.3.1, körtider i sekunder.
körning 1
körning 2
körning 3
medeltal
script1
38,214
38,219
38,277
38,237
script2
36,458
36,434
36,477
36,456
script2b
23,295
23,363
23,371
23,343
I tabell 6.2.2 presenteras körtiderna efter att optimeringarna enligt kapitel
5.3.3 tillämpats. Jämfört med originalversionen har nu medeltalet av
körtiderna för script1 minskat med en faktor på ca. 1,40. För script2 och
script2b
var
inte
optimeringsfaktorn
lika
hög
eftersom
hoppelimineringsmetoden inte täckte hela funktionen algn_fill_row_aff. En
körning av script2 gav en förbättringsfaktor på ca. 1,13 medan motsvarande
förbättring av script2b endast gav en faktor på ca. 1,08.
67
Tabell 6.2.4. Optimering enligt kapitel 5.3.4, körtider i sekunder
körning 1
körning 2
körning 3
medeltal
script1
34,890
34,767
34,739
34,799
I tabell 6.2.4 presenteras tiderna för körningar av POY användande den
slutgiltigt
optimerade
funktionen
algn_fill_row.
Eftersom
endast
algn_fill_row var föremål för denna optimering är endast körtiderna för
script1 av intresse. En 8-gångers upprullning av funktionens slinga minskar
medelkörtiden med ytterligare 4 sekunder jämfört med resultat från
föregående optimering. Den slutgiltiga förbättringsfaktorn blir då ca. 1,54.
Tabell 6.2.5. Optimeringar enligt kapitel 5.3.5 , körtider i sekunder
körning 1
körning 2
körning 3
medeltal
script2
32,367
32,361
32,337
32,355
script2b
21,809
21,799
21,804
21,804
I tabell 6.2.5 presenteras resultaten av de sista optimeringarna (se kapitel
5.3.3.5) för funktionen algn_fill_row_aff. För exekvering av programmet
körande script2 har en uppsnabbningsfaktor på 1,27 erhållits jämfört med
originalversionen. För script2b ligger motsvarande faktor på ca. 1,16.
Figur 6.2.1 samt figur 6.2.1 är profileringar av POY baserade på de mest optimala
versionerna av algn_fill_row och algn_fill_row_aff. Motsvarande profileringar för de
icke optimerade versionerna av funktionerna behandlades i kapitel 5.2. Enligt
resultaten producerade av gprof, har exekveringstiden för algn_fill_row minskat från
72,12 % till 54,96 % och för algn_fill_row_aff har motsvarande tid minskat från 77
% till 69,96 %. Enligt gprof har tiden spenderad i algn_fill_row sjunkit från 34,36s till
17,08s, vilket betyder att algn_fill_row har snabbats upp med en faktor på ca 2,0
(34,36 / 17,08 = 2,01). Motsvarande för algn_fill_row_aff, har tiden den för
funktionen spenderade tiden minskat från 28,12 till 20,51, vilket ger en
uppsnabbningsfaktor på ca 1,37 (28,12 / 20,51 = 1,3710).
68
Figur 6.2.1. Profilering av script1, POY optimerad.
Figur 6.2.2. Profilering av script2, POY optimerad.
Optimering av funktionen algn_fill_row gav bättre resultat än optimering av
algn_fill_row_aff pga. att bättre avgränsade optimeringstekniker kunde
användas. Som exempel kunde nästan alla villkorliga hopp elimineras i den
förstnämnda
funktionen.
Även
om
samma
metod
applicerades
på
algn_fill_row_aff, kunde bara en del av funktionen täckas, varför den relativa
optimeringen blev lägre. Metoden att rulla upp slingor kunde inte heller
tillämpas i algn_fill_row_aff eftersom koden var för komplex för att effektivt
kunna rullas ut.
6.3.
Resultat baserade på mätningar gjorda med PAPI
PAPI (Performance Application Programming Interface, se [DLMMT]) är ett
gränssnitt för mätningar av antal händelser (events) som utförts för olika delar
av
hårdvaran
i
processorn.
Nyare
processorer
är
utrustade
med
hårdvaruräknare (hardware counters) och t.ex. Core 2 Duo-processorn, som
69
nyttjats för de presenterade mätningarna, har två stycken räknare per kärna.
Med hjälp av PAPI kan man specificera händelser för hårdvaruräknarna att
räkna. En händelse kan t.ex. vara antalet L1-cachemissar som uppstår eller
antal villkorliga hopp som utförs under en körning av ett program.
Med PAPI har mätningar gjorts för de mest optimerade versionerna av
algn_fill_row och algn_fill_row_aff samt för motsvarande originalversioner.
Mätningarna har gjorts genom att starta hårdvaruräknarna före anrop till
funktionen i fråga. När funktionen exekverat och återvänder till den funktion
som anropade den, pausas hårdvaruräknarna. På detta sätt kan man avgränsa
sig till att endast mäta händelser i funktionerna. Resultaten nedan baserar sig
på fem miljoner funktionsanrop till respektive funktion.
Tabell 6.3.1. PAPI-mätningar för algn_fill_row.
L2 cachemissar
L2 cacheåtkomster
kvot
original
3 800
12 400 000
0,0003
optimerad
3 800
11 600 000
0,0003
L1 cachemissar
L1 cacheåtkomster
kvot
5 602 730
1 720 000 000
0,0033
5 540 000
1 470 000 000
0,0038
127 000 000
901 000 000
0,1410
20 100 000
205 000 000
0,0980
Fel förutsagda hopp
villkorliga hopp
kvot
I tabell 6.2.4 presenteras olika mätningar av händelser för algn_fill_row. För
den icke optimerade versionen är antalet tillträden till L2-cacheminnet dryga
12 miljoner, av vilka endast 3800 accesser motsvarar cachemissar. Detta
betyder att endast 0,3 % av L2-tillträdena är missar. Om antalet L2-missar
jämförs med t.ex. antalet tillträden till L1-cacheminnet, kan man dra den
slutsatsen att L2-cachemissarna så gott som inte påverkar exekveringstiden av
programmet.
Som
tidigare
nämndes
i
kapitel
5.3.2,
fungerar
minneshanteringen relativt bra i POY. Vidare kan man se att det sker totalt
mindre cachetillträden i den optimerade versionen. Detta kan möjligtvis bero
på att data hålls bättre i registren efter optimeringen. Ett problem som
observerats vid mätningarna är att antalet cachetillträden för L2 är ungefär
70
dubbelt fler än antalet L1-missar. Eftersom en L1-cachemiss, enligt
principerna för cacheminnen i en Core 2 Duo-processor, resulterar i en L2cacheaccess borde antalen motsvara varandra. Vad detta beror på har inte retts
ut.
Den märkbaraste optimeringen av algn_fill_row är minskningen av antalet
villkorliga hopp. I den optimerade versionen enligt tabell 6.2.4 sker drygt fyra
gånger mindre villkorliga hopp och procentuellt har även antalet fel
förutsagda hopp minskat. Orsaken till denna minskning torde vara optimering
5.3.3.1.
Tabell 6.3.2. PAPI-mätningar för algn_fill_row.
L2 cachemissar
L2 cacheåtkomster
kvot
original
7 800
14 100 000
0,0006
optimerad
7 400
11 600 000
0,0006
L1 cachemissar
L1 cacheåtkomster
kvot
10 502 807
5 460 000 000
0,0019
8 620 000
3 950 000 000
0,0022
Fel förutsagda hopp
villkorliga hopp
kvot
110 000 000
2 100 000 000
0,0524
49 500 000
924 000 000
0,0536
I tabell 6.3.2 presenteras de mätningar som gjorts för algn_fill_row_aff. Likt
resultaten
i
tabell
6.3.1
är
L2-cachemissarna obetydliga.
Även
L1-
cachemissarna är så få att de knappast nämnvärt påverkar körtiden. Antalet
L1-cacheminnesåtkomster har däremot minskat från ca 5,5 miljarder till
knappt 4,0 miljarder. Detta beror antagligen främst på de optimeringar som
tillämpades i Kapitel 5.3.5. Även antalet villkorliga hopp har halverats, vilket
främst torde vara ett resultat av optimering 5.3.3.1.
71
7. Avslutning
Inom fylogenin strävar man efter att gruppera organismer och ordna dessa
enligt relativ släktskap. Kladogram, som är en form av specialiserade binära
träd, används ofta för visualisering av hypotetisk släktskap organismerna
sinsemellan. Som data för dessa hypotetiska släktskapsbildningar fungerar ur
organismerna erhållna DNA-sekvenser. Dessa sekvenser har utvunnits i
laboratoriemiljöer och avkodats till textsträngar som datorer kan behandla.
Med datorns hjälp, i form av program som bearbetar dylika data, har man
kunnat återskapa verklighetstrogna evolutionsförlopp.
POY är ett program som utnyttjar en teknik kallad dynamisk homologi
(Dynamic Homology) för skapandet av kladogram. Huvudprincipen bakom
denna teknik grundar sig i en strävan att välja det bästa av alla unika möjliga
kladogram
som
teoretiskt
kan
bildas
med
antalet
organismer
som
utgångspunkt. För varje hypotetiskt kladogram som granskas, kan olika
metoder appliceras för testning av dess verklighetstrogenhet. En av de
viktigaste metoderna kallas för direkt optimering (Direct Optimization) och är
den för tillfället huvudsakligen mest använda metoden.
Att med hjälp av POY bilda stora men även pålitliga kladogram, baserade på
hundratals till tusentals organismer med DNA-sekvenser på tusentals element,
kan ta veckor eller månader att exekvera beroende på hur noga sökningarna
är. Att vänta så långa tider kan för användaren vara frustrerande, varför en
optimering av programmet är relevant. Att optimera ett program innebär att
det görs snabbare med hjälp av någon, ofta känd, optimeringsteknik.
Optimering på kodnivå, som detta arbete främst har handlat om, grundar sig i
hur man ska gå tillväga då man skriver kod för att kunna utnyttja
processorkraften maximalt. I avhandlingen har läsaren även givits en inblick i
hur moderna processorer fungerar. Att känna till de viktigaste egenskaperna
hos moderna processorer är en väsentlig del vad gäller kodoptimering.
72
Målet med arbetet som gjorts har varit att försöka snabba upp programmet så
mycket som möjligt, vilket nåtts till en viss del. I denna avhandling har
läsaren bekantat sig med olika kodoptimeringsmetoder som t.ex. upprullning
av slingor (kapitel 5.3.4) och eliminering av villkorliga hopp (kapitel 5.3.3).
Dessa metoder har applicerats i försök att göra olika delar av POY snabbare. I
kapitel sex behandlas resultaten av optimeringarna. Det visade sig att
programmet var tämligen bra optimerat redan från början, varför tillämpning
av optimeringsmetoderna presenterade i denna avhandling endast har
inbringat en maximal optimeringsfaktor på drygt 1,5.
73
Källor
[AMD04] Advanced Micro Devices, Inc. (2004). Software Optimization Guide for
AMD Athlon™ 64 and AMD Opteron™ Processors.
[Bray57] Brayer E.F. (1957), Calculating the Standard Error of a Proportion.
Applied Statistics, (6)1: 67 – 68.
[Casa07] Casanova, H. (2007). Lecture Notes. The x86 Architecture. Machine-Level
and Systems Programming. Information and Computer Sciences Department,
University of Hawaii.
navet.ics.hawaii.edu/~casanova/courses/ics312_fall07/slides/ics312_x86.pdf
[Cpui06] (2006). 486. En kort artikel från WWW.CPU-INFO.COM om x486processorn. www.cpu-info.com/index2.php?mainid=html/cpu/486.php
[Cole] Coleman, C.L. (Årtal okänt). En plagierad sammansättning. Using Inline
Assembly With gcc. Free Software Foundation, Inc., Boston, USA.
www.cs.virginia.edu/~clc5q/gcc-inline-asm.pdf
[DLMMT] Dongarra, J., London, K., Moore, S., Mucci, P., Terpstra, D. (2001). Using
PAPI for hardware performance monitoring on Linux systems. Conference on Linux
Clusters: The HPC Revolution, Linux Clusters Institute, Urbana, Illinois
[FaHa03] Farber, D.A., Hall, B. (2003). The problematics of the Pareto Principle.
Boalt Working Papers in Public Law (35), University of California, Berkley.
[Farr70] Farris, J.S. (1970). Methods for computing Wagner trees. Systematic
Zoology 19: 83 – 92.
[Fels78] Felsenstein, J. (1978). The number of evolutionary trees. Systematic
Zoology, 27: 27 – 33.
[FeSt00] Fenlason, J. & Stallman, R. (2000). GNU gprof, The GNU Profiler. Free
Software Foundation, Inc., Boston, USA. www.skyfree.org/linux/references/gprof.pdf
[Fog07a] Fog, A. (2007). Optimizing software in C++ – An optimization guide for
Windows, Linux and Max platforms. Copenhagen University College of Engineering.
[Fog07b] Fog, A. (2007). The microarchitecture of Intel and AMD CPU’s – An
optimization guide for assembly programmers and compiler makers. Copenhagen
University College of Engineering.
[Gatl01] Gatliff, B. (2001). Getting Started with GNU:A Tutorial Introduction Using
the ARM Evaluator-7T. www.microcross.com/gnu-arm7t-microcross.
[Glas04] Glaskowsky, P.N. (2004). Prescott Pushes Pipelining Limits.
MICROPROCESSOR REPORT. Reed Electronics Group.
74
[Hand07] Okänd skribent. (2007). A Handbook, Help Me Understand Genetics.
Lister Hill National Center for Biomedical Communications U.S. National Library of
Medicine National Institutes of Health. ghr.nlm.nih.gov/handbook.pdf
[Inte06] Intel Corp. (2006). Intel® 64 and IA-32 Architectures Software Developer’s
Manual, vol.1 Basic Architecture.
[KYB03] Korf, I., Yandell, M., Bedell, J. (2003). ”BLAST”: An Essential Guide to
the Basic Local Alignment Search Tool.
[LeSm84] Lee, J.K.L., Smith, A.J. (1984). Branch prediction strategies and branch
target buffer design. IEEE Computer, 17(1), 6 - 22.
[LeWi00] Leeuwen, J. & Wiedermann J., (2000.) The Turing Machine Paradigm in
Contemporary Computing.Technical Report. Department of Computer Science,
Utrecht University, Institute of Computer Science, Academy of Sciences of the Czech
Republic.
[Loh05] Loh, G. (2005). Memory and Caches. Lecture notes. Georgia Institute of
Technology. www-static.cc.gatech.edu/classes/AY2005/cs3220_spring/cachenotes.pdf
[Mnei02a] Mneimneh, S. (2002). Lecture Notes. Lecture 5: Affine gap penalty
function, multiple sequence alignment. Massachusetts Institute of Technology.
engr.smu.edu/~saad/courses/cse8354/lectures/lecture5.pdf
[Mnei02b] Mneimneh, S. (2002). Lecture Notes. Lecture 6: Affine gap penalty
function, multiple sequence alignment. Massachusetts Institute of Technology.
engr.smu.edu/~saad/courses/cse8354/lectures/lecture6.pdf
[Moor65] Moore, G.E. (1965). Cramming more components onto integrated circuits.
Electronics, 38(8).
[NeWu70] Needleman, S., Wunsch, C. (1970). A general method applicable to the
search for similarities in the amino acid sequence of two proteins. J. Mol. Biol.,
48(3): 443 – 453.
[NoCa77] Platnick, N. & Cameron, H., (1977). Cladistic methods in textual,
linguistic, and phylogenetic analysis. Systematic Zoology, 26: 380 – 385.
[Pas02] Pas, R. (2002). Memory Hierarchy in Cache-Based Systems. High
Performance Computing, Sun Microsystems, inc., Network Circle, Santa Clara,
California, USA.
[Patt03] Patterson, J. (2003). Modern Microprocessors A 90 Minute Guide!
www.pattosoft.com.au/Articles/ModernMicroprocessors/
[Pime] Neeman, H. (Årtal okänt). A syllabus to Introduction to Parallel Architecture.
(undervisningsmaterial). Computer Systems Architecture group, Informatics Institute,
University of Amsterdam. http://staff.science.uva.nl/~andy/aci/syl.pdf
75
[PrEn97] Preiss, B.R., Eng, P. (1997). Data Structures and Algorithms with ObjectOriented Design Patterns inC++, Complete Trees.
www.brpreiss.com/books/opus4/html/page355.html
[PTVF07] Press, W., Teukolsky, S., Vetterling, T., Flannery, B. (2007). Numerical
Recipes 3rd Edition: The Art of Scientific Computing.
[RPF03] Rico, R., Pérez, J-I., Frutos A. (2003). The impact of x86 instruction set
architecture on superscalar processing. Department of Computer Engineering,
Universidad de Alcala´ , Alcala de Henares, Spain.
[Smit82] Smith, A.J. (1982). Cache Memories. Computing Surveys 14(3). University
of California, Berkley.
[Stal03] Stallman, R. & Co. (2003). Using the GNU Compiler Collection (GCC). Free
Software Foundation, Inc. gcc.gnu.org/onlinedocs/gcc-4.1.2/gcc.pdf
[Tola01] Toland, E. (2001)., DNA Mutations. Genetics 101.
www.genetichealth.com/G101_Changes_in_DNA.shtml
[Torr07] Torres, G. (2007). How The Memory Cache Works. Hardware Secrets Web
Tutorial. www.hardwaresecrets.com/article/481/
[VVBW] Varón, A., Vinh, L.S., Bomash, I., Wheeler, W.C. (2007). Användarmanual,
POY 4.0 Beta, Quickstart Guide and Program Documentation. American Museum of
Natural History, New York.
[Whee96] Wheeler, W. (1996). Optimization alignment: the end of multiple sequence
alignment in phylogenetics?. Department of Invertebrates, American Museum of
Natural History.
[Whee99] Wheeler, W.C. (1999). Fixed character states and the optimization of
molecular sequence data. Cladistics 15: 379 – 385.
[Whee01] Wheeler, W.C. (2001). Homology and the optimization of DNA sequence
data. Cladistics 17: 3 – 11.
[Whee03a] Wheeler, W.C. (2003). Iterative pass optimization of sequence data.
Cladistics 19: 254–260.
[Whee03b] Wheeler, W.C. (2003). Search-based optimization. Cladistics 19: 348 –
355.
[Whee06] Wheeler, W.C (2006). Dynamic Homology and Phylogenetic Systematics:
A Unified Approach Using POY. American Museum of Natural History, New York.
[WSBF91] Whiley, E.O., Siegel-Causey, D., Brooks D.R., Funk V.A. (1991) The
Compleat Cladist A Primer of Phylogenetic Procedures. Special publication No. 19.
University of Kansas Museum of Natural History.
76
77