Skip to content

Encode horizontal deltas using +1/-1 indicators for up to 20% speed gain #214

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
RagnarGrootKoerkamp wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Encode horizontal deltas using +1/-1 indicators for up to 20% speed gain #214

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 31 additions & 18 deletions edlib/src/edlib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -401,46 +401,45 @@ static inline unsigned char* createReverseCopy(const unsigned char* const seq, c
* @param [in] Pv Bitset, Pv[i] == 1 if vin is +1, otherwise Pv[i] == 0.
* @param [in] Mv Bitset, Mv[i] == 1 if vin is -1, otherwise Mv[i] == 0.
* @param [in] Eq Bitset, Eq[i] == 1 if match, 0 if mismatch.
* @param [in] hin Will be +1, 0 or -1.
* @param [in] Phin Bitset, 00..01 when hin == +1.
* @param [in] Mhin Bitset, 00..01 when hin == -1.
Comment on lines +404 to +405

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* @param [in] Phin Bitset, 00..01 when hin == +1.
* @param [in] Mhin Bitset, 00..01 when hin == -1.
* @param [in] Phin Bitset, 00..01 when hin == +1, otherwise 00..00.
* @param [in] Mhin Bitset, 00..01 when hin == -1, otherwise 00..00.

* @param [out] PvOut Bitset, PvOut[i] == 1 if vout is +1, otherwise PvOut[i] == 0.
* @param [out] MvOut Bitset, MvOut[i] == 1 if vout is -1, otherwise MvOut[i] == 0.
* @param [out] hout Will be +1, 0 or -1.
* @param [out] Phout Bitset, 00..01 when hout == +1.
* @param [out] Mhout Bitset, 00..01 when hout == -1.
Comment on lines +408 to +409

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is actually a single output, a pair, not two outputs, so let's show it as such here, as a single @param?

Comment on lines +408 to +409

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* @param [out] Phout Bitset, 00..01 when hout == +1.
* @param [out] Mhout Bitset, 00..01 when hout == -1.
* @param [out] Phout Bitset, 00..01 when hout == +1, otherwise 00..00.
* @param [out] Mhout Bitset, 00..01 when hout == -1, otherwise 00..00.

*/
static inline int calculateBlock(Word Pv, Word Mv, Word Eq, const int hin,
static inline pair<Word, Word> calculateBlock(Word Pv, Word Mv, Word Eq, const Word Phin, const Word Mhin,
Word &PvOut, Word &MvOut) {
// hin can be 1, -1 or 0.
// 1 -> 00...01
// 0 -> 00...00
// -1 -> 11...11 (2-complement)
Comment on lines 413 to 416

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not that meaningful any more, we can maybe rewrite it like this:

Suggested change
// hin can be 1, -1 or 0.
// 1 -> 00...01
// 0 -> 00...00
// -1 -> 11...11 (2-complement)
// (Phin, Mhin) can have one of three possible values:
// - (1, 0) if hin is 1.
// - (0, 1) if hin is -1.
// - (0, 0) if hin is 0.


Word hinIsNeg = static_cast<Word>(hin >> 2) & WORD_1; // 00...001 if hin is -1, 00...000 if 0 or 1

Word Xv = Eq | Mv;
// This is instruction below written using 'if': if (hin < 0) Eq |= (Word)1;
Eq |= hinIsNeg;
Eq |= Mhin;
Word Xh = (((Eq & Pv) + Pv) ^ Pv) | Eq;

Word Ph = Mv | ~(Xh | Pv);
Word Mh = Pv & Xh;

int hout = 0;
// This is instruction below written using 'if': if (Ph & HIGH_BIT_MASK) hout = 1;
hout = (Ph & HIGH_BIT_MASK) >> (WORD_SIZE - 1);
Word Phout = Ph >> (WORD_SIZE - 1);

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see that you removed & HIGH_BIT_MASK, and I was trying to figure why, but now that I am looking at it, it seems it was never useful, was it? So this was always redundant not a change caused by the other changes here, right?
Same for Mhout below.

// This is instruction below written using 'if': if (Mh & HIGH_BIT_MASK) hout = -1;
hout -= (Mh & HIGH_BIT_MASK) >> (WORD_SIZE - 1);
Word Mhout = Mh >> (WORD_SIZE - 1);

Ph <<= 1;
Mh <<= 1;

// This is instruction below written using 'if': if (hin < 0) Mh |= (Word)1;
Mh |= hinIsNeg;
Mh |= Mhin;
// This is instruction below written using 'if': if (hin > 0) Ph |= (Word)1;
Ph |= static_cast<Word>((hin + 1) >> 1);
Ph |= Phin;

PvOut = Mh | ~(Xv | Ph);
MvOut = Ph & Xv;

return hout;
return {Phout, Mhout};
}

/**
Expand Down Expand Up @@ -581,17 +580,22 @@ static int myersCalcEditDistanceSemiGlobal(

int bestScore = -1;
vector<int> positions; // TODO: Maybe put this on heap?
const int startHout = mode == EDLIB_MODE_HW ? 0 : 1; // If 0 then gap before query is not penalized;
const Word PstartHout = mode == EDLIB_MODE_HW ? 0 : 1; // If 0 then gap before query is not penalized;

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think now is a bit harder to track what is start state, since only P part of hout is set at start, and what it the relationship between the start hout and Phout and Mhout.

What if we did this instead:

    int bestScore = -1;
    vector<int> positions; // TODO: Maybe put this on heap?
    // PMhout is (Phout, Mhout): (1,0) if hout is 1, (0,1) if hout is -1, (0,0) if hout is 0.
    // Also, hout == Phout - Mhout.
    const pair<Word,Word> startPMhout = mode == EDLIB_MODE_HW ? {0,0} : {1,0}; // If 0 then gap before query is not penalized;
    const unsigned char* targetChar = target;
    for (int c = 0; c < targetLength; c++) { // for each column
        const Word* Peq_c = Peq + (*targetChar) * maxNumBlocks;

        //----------------------- Calculate column -------------------------//
        pair<Word,Word> PMhout = startPMhout;
        bl = blocks + firstBlock;
        Peq_c += firstBlock;
        for (int b = firstBlock; b <= lastBlock; b++) {
            PMhout = calculateBlock(bl->P, bl->M, *Peq_c, PMhout.first, PMhout.second, bl->P, bl->M);
            bl->score += PMhout.first - PMhout.second; // PMhout.first - PMhout.second == hout
            bl++; Peq_c++;
        }
        bl--; Peq_c--;
        const int hout = PMhout.first - PMhout.second;
        //------------------------------------------------------------------//

So now we have startPMhout at the start, as a single value that describes starting condition, then we have PMhout to track how hout develops as we calculate block by block, and I avoided recalculating hout in every step, we don't really need it in that format, PMhout is actually the format we need. Only at the very end we actually calculate hout, because it is useful for the readability further below.

const unsigned char* targetChar = target;
for (int c = 0; c < targetLength; c++) { // for each column
const Word* Peq_c = Peq + (*targetChar) * maxNumBlocks;

//----------------------- Calculate column -------------------------//
int hout = startHout;
Word Phout = PstartHout;
Word Mhout = 0;
int hout = 0;
bl = blocks + firstBlock;
Peq_c += firstBlock;
for (int b = firstBlock; b <= lastBlock; b++) {
hout = calculateBlock(bl->P, bl->M, *Peq_c, hout, bl->P, bl->M);
pair<Word,Word> PMhout = calculateBlock(bl->P, bl->M, *Peq_c, Phout, Mhout, bl->P, bl->M);
Phout = PMhout.first;
Mhout = PMhout.second;
hout = Phout - Mhout;
bl->score += hout;
bl++; Peq_c++;
}
Expand All @@ -605,7 +609,8 @@ static int myersCalcEditDistanceSemiGlobal(
lastBlock++; bl++; Peq_c++;
bl->P = static_cast<Word>(-1); // All 1s
bl->M = static_cast<Word>(0);
bl->score = (bl - 1)->score - hout + WORD_SIZE + calculateBlock(bl->P, bl->M, *Peq_c, hout, bl->P, bl->M);
pair<Word, Word> PMhout = calculateBlock(bl->P, bl->M, *Peq_c, Phout, Mhout, bl->P, bl->M);
bl->score = (bl - 1)->score - hout + WORD_SIZE + PMhout.first - PMhout.second;
Comment on lines +612 to +613

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we apply the changes I suggested above, then we also need to adjust this part a bit:

Suggested change
pair<Word, Word> PMhout = calculateBlock(bl->P, bl->M, *Peq_c, Phout, Mhout, bl->P, bl->M);
bl->score = (bl - 1)->score - hout + WORD_SIZE + PMhout.first - PMhout.second;
pair<Word, Word> nextPMhout = calculateBlock(bl->P, bl->M, *Peq_c, PMhout.first, PMhout.second, bl->P, bl->M);
const int nextHout = nextPMhout.first - nextPMhout.second;
bl->score = (bl - 1)->score - hout + WORD_SIZE + nextHout;

} else {
while (lastBlock >= firstBlock && bl->score >= k + WORD_SIZE) {
lastBlock--; bl--; Peq_c--;
Expand Down Expand Up @@ -780,9 +785,14 @@ static int myersCalcEditDistanceNW(const Word* const Peq, const int W, const int

//----------------------- Calculate column -------------------------//
int hout = 1;
Word Phout = 1;

@Martinsos Martinsos May 15, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If changes I suggested above make sense, we could do the same thing here -> have only PMhout here, no Phout, Mhout, and hout, till we are done with calculating blocks, only then we calculate hout, and below we go with nextPMhout and nextHout.
Oh I see now I used newHout here already -> but let's rename that to nextHout, that is a better name (since it is "next" block).

Word Mhout = 0;
bl = blocks + firstBlock;
for (int b = firstBlock; b <= lastBlock; b++) {
hout = calculateBlock(bl->P, bl->M, Peq_c[b], hout, bl->P, bl->M);
pair<Word, Word> PMhout = calculateBlock(bl->P, bl->M, Peq_c[b], Phout, Mhout, bl->P, bl->M);
Phout = PMhout.first;
Mhout = PMhout.second;
hout = Phout - Mhout;
bl->score += hout;
bl++;
}
Expand All @@ -806,7 +816,10 @@ static int myersCalcEditDistanceNW(const Word* const Peq, const int W, const int
lastBlock++; bl++;
bl->P = static_cast<Word>(-1); // All 1s
bl->M = static_cast<Word>(0);
int newHout = calculateBlock(bl->P, bl->M, Peq_c[lastBlock], hout, bl->P, bl->M);
pair<Word, Word> PMhout = calculateBlock(bl->P, bl->M, Peq_c[lastBlock], Phout, Mhout, bl->P, bl->M);
Phout = PMhout.first;
Mhout = PMhout.second;
int newHout = Phout - Mhout;
bl->score = (bl - 1)->score - hout + WORD_SIZE + newHout;
hout = newHout;
}
Expand Down