Very cool! Tell me when you need a code review

Agree that the first version should be naive. Can improve from there!

Hey @TwFlem thanks for submitting this draft!

For naive I meant really, really naive. I'm not sure what everything here is doing yet but here's a simple 30 line BWT implementation that works:

func burrowsWheelerTransform(input string) string {
	// Add a sentinel character to the end of the input
	input += "$"

	// Create a matrix to store rotations
	rotations := make([]string, len(input))
	for i := 0; i < len(input); i++ {
		rotations[i] = input[i:] + input[:i]
	}

	// Sort the rotations lexicographically
	sort.Strings(rotations)

	// Extract the last characters of each rotation to get the transformed string
	var result string
	for _, rotation := range rotations {
		result += string(rotation[len(rotation)-1])
	}

	return result
}

Obviously yours is more complex and very likely faster but I'm not sure what's going on yet. What's the best way to compare this^^ toy version with your own in this PR? What references did you use to create your draft? Any links?

-Tim

Sure! So I read through many materials, but the most salient of them all came from Ben Langmead. He has a whole youtube playlist about the subject.

I wish it wasn't so complex! So the BWT by itself is actually a small piece of the puzzle- there are auxiliary components that go along with the BWT to make it possible to both compress the size from the original string while also making it usable as a search index to perform other kinds of analysis like sequence alignment.

The equivalent to what you wrote above is here in the pr (ignoring the skiplist and the part about wavelet tree). The only difference with what I have here is that it saves some memory because we don't have to store NxN matrix of the string before we sort it. The getBWTIndex is what allows us to get that original character from the provided sequence that would show up in the L column for that given rotation.

Everything else in this PR is what is necessary to enable search and alignment.

To sum up what I have in the PR and why it is there at a high level:

1. Bitvector- gives us a way to encode information in a memory efficient way. For example, instead of needing 1 byte to represent a character, all we need it 1 bit in the bitvector.
2. RSABitvector- RSA stands for Rank, Select, Access. So far, only Rank and Access are needed in the current implementation. However, select is needed for the RLFM/r-index implementations which would be the most usable implementation for anyone trying to do any serious work with the BWT. More on RLFM and r-index later though.
3. Wavelet tree- This is a data structure that allows us to make RSA queries on strings. TL;DR it is a binary tree that allows us to sub divide a given string and corresponding alphabet so we can make RSA queries. This is probably the most complex thing. I suggest just giving this a quick once over.
4. LF mapping: LF stands for Last and First: Last column of the BWT transform and First column of the BWT transform. In this PR, F is the skiplist and L is just l for now. This is actually how we conduct search and alignment. In the Count, Locate, and Extract queries, you'll see we start with the whole range of both coulmns as the search range and use the skiplist along with ranking the L column to narrow down that search range. I suggest watching this to see how the backwards search is conducted. It's a lot easier to see what's going on in the visual.
5. prefixArray: This is a naive implementation. This takes up way too much memory as it currently is. The only reason why I didn't improve it further is because we'd most likely have to completely change its implementation if we want to do the RLFM and or r-index. The prefix array gives us a look up table that allows us map characterInTheFCoulmn -> characterInTheOriginalString. This allows us to extract the text.

Also, I mentioned the RLFM and r-index. The only reason why I didn't continue with doing these is because this has gotten complex enough... that and I figured it'd be good to make sure you want all of this in the first place. With that being said, I highly recommend we pursue them. In a very similar manner to how the skipList in this PR is constructed, we can compress the L column substantially- especially for sequences that are very repetitive. He has some very compelling numbers in this video towards the end. Also, once you see the compression, it makes sense how much of an improvement this could be.

Btw, I will fix the names and do some refactoring before taking this out of draft :).

It'd be cool to do a little talk and demo for this.

I'd be happy to do that. Right now, there's not too much to demo or talk about maybe. We can talk about saving memory, indexing strings, the magic of the BWT and how rank allows us to take advantage of that magic.

Maybe it'd be even more valuable to have a demo once we have the fuzzy matching that Koeng101 mentioned in #360 along with having the RLFM and r-index? At that point, the BWT seems like it would be feature complete enough to reach a critical mass of usefulness for Poly's users.

It's a great example for teaching and maybe we should share it back with Ben? Just some late night thoughts.

I'd definitely like to at least thank him! It's amazing that he is giving that knowledge to the community like that. Perhaps his channel deserves a spot in https://github.com/TimothyStiles/how-to-synbio. I'm happy to share anything else that he might appreciate as well!

It's not clear where a lot of the test cases are coming from and there should be more comments about where they're from, and what corner case their testing.

I added some comments about some edge cases that are being tested like with Jacobson's rank and the rsa_bitvector rank tests. I added a mention of how useful the the reconstruction test is for the WaveletTree. I also added a reconstruction test that pretty much guarantees that whatever we use for the suffix array in the future is well formed- otherwise we would never be able to reconstruct the string.

I added some more tests and more validation. So now we throw errors for empty and other invalid inputs. I fixed an edge case where a user could supply a string with a single alphabet character.

Otherwise, there aren't really too many edge cases to go around with the FM index BWT. It's all predicated on the fact that the LF search works. LF search is all predicated on the fact that the Wavelet Tree and other auxiliary DSAs work as expected. The other tests are just to meant the show that auxiliary pieces work as expected. The RSA and Bit vector test cases are verbose because of all the bit math involved. Those are there to give us peace of mind that it works as expected.

There will be some more BWT edge cases with the other improvements that I'll follow up with like the RLFM index (done) and the r-index (in progress). Those I can definitely call out. Those changes will be pretty small btw. Nothing like getting the basic BWT off the ground 😭

I'm with @carreter that it'd be good to wrap my head around the actual concepts but I'm also thinking we should find some outside folk more familiar with bwt to review this too if we can.

Sounds good! It might be better to emphasize the correctness over performance improvements. Reasonable performance improvements are already on the way 👀