# Lucas van Dijk > Personal blog on high-performance software engineering to solve biological problems at scale. Public Ghost content for AI and LLM tooling. This file includes a bounded export of public pages first, then recent public posts. Append `.md` to any post or page URL to get the content in Markdown (for example, `/example-post.md`). ## Pages ### About URL: https://lucasvandijk.nl/about/ Last updated: 2025-11-23T23:59:37.000Z I am a scientist and research software engineer with 7+ years of experience, passionate about developing high-performance tools and machine learning systems to accelerate research. My tools and algorithms have enabled new insights into human disease and enabled genome analyses previously impossible. ## Areas of Expertise - **Languages & Frameworks.** Python, C++, Rust, JavaScript, scikit-learn, TensorFlow, Jax, Pandas, NumPy - **Bioinformatics.** Sequence Alignment, Variant Calling, Pangenomics, Genome Assembly, Phylogenetics. - **Machine Learning.** Deep Learning, Hidden Markov Models, Logistic Regression, Linear Models - **Software Engineering.** Git, Docker, CI/CD, Unit Testing, Google Cloud Platform, PostgreSQL, Linux ## Professional experience ***Pacific Biosciences, Remote*** - Senior Software Engineer, Bioinformatics, Nov 2025 - Current ***Broad Institute of MIT and Harvard, Cambridge, MA*** - Computational Scientist II, *Data Sciences Platform,* Sep 2024 - Jul 2025 - Developed a Google Cloud-based data pre-processing pipeline, curating disease phenotypes from electronic health records (EHR) for 10,000 participants in the “All of Us” biobank and assessing the quality of associated terabyte-scale next-generation sequencing (NGS) data. - Led the design of interpretable machine learning models integrating genomic (PacBio HiFi, Illumina), transcriptomic, and proteomic data to predict disease onset and identify novel biomarkers. - Redesigned and optimized an implementation of a hidden markov model for recombination-aware DNA sequence alignment using Cython and C++, resulting in 360x faster inference. [GitHub - castcollab/tesserae2: Tesserae2: Fast recombination-aware global and local alignment.Tesserae2: Fast recombination-aware global and local alignment. - castcollab/tesserae2![](https://lucasvandijk.nl/content/images/icon/pinned-octocat-093da3e6fa40.svg)GitHubcastcollab![](https://lucasvandijk.nl/content/images/thumbnail/tesserae2)](https://github.com/castcollab/tesserae2?ref=lucasvandijk.nl) - Computational Associate II, *Bacterial Genomics Lab,* Nov 2017 - Aug 2024 - Designed and implemented a new optimal partial order alignment algorithm in Rust, which was, on average, 4.1x faster than similar tools, enabling DNA sequence-to-graph alignments previously impossible. - Led the development of Python and C++-based software specifically designed to track low-abundance (>0.1%) bacterial strains in complex microbial communities (e.g., the human gut microbiome) using *whole metagenome sequencing* data. - Developed a cloud-based pipeline using the Workflow Description Language (WDL) and Docker to enable fast and reproducible characterization of thousands of metagenomic samples. - Obtained detailed insight into the *E. coli* strain-level dynamics in a year-long longitudinal microbiome study of women with recurrent urinary tract infections (UTIs), revealing unexpected similarities with a healthy control group and that the UTI-causing strain is rarely cleared from the gut after antibiotics. [StrainGE: Strain Genome ExplorerStrainGE is a toolkit for tracking and characterizing low-abundance strains in complex microbial communities. It enables detailed insights into the bacterial strain-level diversity of whole metagenomic sequencing samples.![](https://lucasvandijk.nl/content/images/icon/favicon-3.ico)Lucas van DijkLucas van Dijk![](https://lucasvandijk.nl/content/images/thumbnail/strainge-cover-2-1.png)](https://lucasvandijk.nl/publications/strainge-strain-genome-explorer/) [Gut microbiome dysbiosis linked with recurrent UTIsMore than half of the women in the US get a urinary tract infection (UTI) in their lifetime, which frequently becomes recurrent. In this paper, we investigated the role of the gut microbiome in facilitating recurrence.![](https://lucasvandijk.nl/content/images/icon/favicon-4.ico)Lucas van DijkLucas van Dijk![](https://lucasvandijk.nl/content/images/thumbnail/umb-nmicro-2.jpg)](https://lucasvandijk.nl/publications/gut-microbiome-dysbiosis-linked-with-recurrent-utis/) ***DSM-firmenich, Delft, The Netherlands*** - Intern, Jun 2016 - Oct 2016 - Reviewed literature and proposed a plan to analyze a large-scale protein production problem with *Bacillus subtilis* using genome-scale metabolic models. ***Studio bereikbaar, Rotterdam, The Netherlands*** - Software Engineer, May 2013 - Jun 2016 - Led the development of web-based geographic information systems (GIS) tools using Django and PostgreSQL, securing multi-million euro infrastructure contracts with the government through improved collaborative project planning and design. - Enabled local analysis and modification of GIS data in the open-source desktop application QGIS by developing a custom Python-based plugin, communicating with a central server through a REST API. ***Thales, Delft, The Netherlands*** - Intern, Sep 2012 - Nov 2012 - Developed a GPU-accelerated tool to analyze the radar reflectivity of navy ships using nVidia's OptiX CUDA ray tracing library. ## Extracurricular Projects ***Python Software Foundation (VisPy Google Summer of Code 2015)*** - Implemented a high-performance graph visualization system in Python and OpenGL, including several automatic graph layout algorithms. - Contributed the open-source code upstream to VisPy, a high-performance scientific data visualization software library. [Drawing arbitrary shapes with OpenGL pointsPart of my Google Summer of Code project involves porting several arrow heads from Glumpy to Vispy. I also want to make a slight change to them: the arrow heads in Glumpy include an arrow body, I want to remove that to make sure you can put an arrow head![](https://lucasvandijk.nl/content/images/icon/favicon-5.ico)Lucas van DijkLucas van Dijk![](https://lucasvandijk.nl/content/images/thumbnail/cover-opengl-1.png)](https://lucasvandijk.nl/2015/06/drawing-arbitrary-shapes-with-opengl-points/) [Home — VisPy![](https://lucasvandijk.nl/content/images/icon/favicon-6.ico)VisPy - Home![](https://lucasvandijk.nl/content/images/thumbnail/galaxy.png)](https://vispy.org/?ref=lucasvandijk.nl) ***Delft University of Technology, "Helios 3D team"*** - Built a large 3x3x1 meter 3D RGB LED cube display as part of a student team. - Architected and implemented the C++ embedded software responsible for receiving the image to display over Wi-Fi and driving individual LEDs to the correct color. ## Education - **PhD Bioinformatics.** Delft University of Technology. 2025. - Research conducted at Broad Institute’s Bacterial Genomics Lab in collaboration with the Delft Bioinformatics Lab. Courses: Immunology, Bayesian Methods for Machine Learning, Deep Learning - **MSc Computer Science.** Delft University of Technology. 2017. - **BEng Electrical Engineering.** The Hague University of Applied Sciences. 2013. ### Sign Up URL: https://lucasvandijk.nl/signup/ Last updated: 2025-03-12T20:37:26.000Z _No content available._ ### Subscribe URL: https://lucasvandijk.nl/subscribe/ Last updated: 2025-03-12T20:37:45.000Z _No content available._ ### Sign In URL: https://lucasvandijk.nl/signin/ Last updated: 2025-03-12T20:38:19.000Z _No content available._ ### Sidebar URL: https://lucasvandijk.nl/sidebar/ Last updated: 2025-11-23T23:48:54.000Z Hi, I'm Lucas! I write about bioinformatics and high-performance software engineering. Expect deep dives into the technical aspects behind a paper or explainers in simple language. Topics I find interesting include sequence alignment algorithms, immunology, and GPU programming. I am a senior software engineer on the Instrument Analysis team at [PacBio](https://pacb.com/?ref=lucasvandijk.nl), a company that builds long-read DNA sequencing machines. We develop high-performance software that transforms sensor data into raw base calls and consensus algorithms that deliver highly accurate DNA reads. Views on this website are my own. [Read the rest of my resume »](https://lucasvandijk.nl/about/) ## Posts ### Relative positional embeddings with RoPE URL: https://lucasvandijk.nl/2025/08/relative-positional-embeddings-with-rope/ Last updated: 2025-08-20T17:36:21.000Z A core component of modern transformer-based language models is self-attention. With self-attention, the embedding of a specific token (or word) in a sentence is transformed by the embeddings of surrounding tokens, capturing sentence context. For example, the word "station" has different meanings in the sentences "I am listening to a radio station" and "I just left the train station". Self-attention by itself is order-agnostic and would compute the same results for the sentence "left train I station the just" or "I just left the train station". This is unfortunate, because word order is important in natural language. So, how do we ensure that self-attention will consider token positions? ## Position as a sinusoidal signal The trick presented in the original "[Attention is all you need](https://arxiv.org/abs/1706.03762?ref=lucasvandijk.nl)" paper is to encode the position of a token into the embedding before computing attention values (Vaswani *et al.* 2017). These position encodings, or *positional embeddings,* can take various forms, and a popular approach is the *sinusoidal* positional embedding. The *d-*dimensional positional encoding for a position *n* is defined as follows: $$ r\_{ni} = \\begin{cases} \\sin(\\frac{n}{L^{i/d}}) & \\text{if } i \\text{ is even}\\\\ \\cos(\\frac{n}{L^{(i-1)/d}}) & \\text{if } i \\text{ is odd}\\end{cases}$$ Let's plot these values for varying *n* and *i* to obtain a better sense of what this encoding represents. ![](https://lucasvandijk.nl/content/images/2025/07/positional-embedding-1.png) The sine (blue) and cosine (orange) waves for different sequence positions (x-axis) and embedding dimensions (y-axis). The encoding values lie in the range \[-1, 1\], with the values in the lower embedding dimensions changing more frequently than those in the higher dimensions. This is similar to binary numbers, where the least significant bits change much more frequently than the most significant bits (Bishop, 2023). This positional embedding vector is then simply added to a token embedding before computing attention: $$ \\textbf{x}\_n' = \\textbf{x}\_n + \\textbf{r}\_n$$ A downside of this approach is that it encodes the *absolute* position into the token embedding. In many sequences, the *relative* ordering of words or tokens is more important. For example, in the sentences "I just left the train station" and "Because I was late, I left the train station," the relationship between "train" and "station" is the same, even though they occur at different absolute positions in each sentence. This is where "Rotary Positional Encodings" (RoPE) come in (Su *et al.* 2023). RoPE similarly uses a sinusoidal signal to encode the position. However, instead of adding this signal to the original input vector, RoPE multiplies the sinusoidal signal with the input after projection into query or key space. For example, to compute a position-adjusted query vector from the original token at position $m$, with the embedding dimension $d=2$: $$ \\textbf{q}' \_m = \\textbf{R} \\textbf{W}\_q \\textbf{x}\_m = \\begin{bmatrix} \\cos m\\theta & -\\sin m\\theta \\\\ \\sin m\\theta & \\cos m\\theta\\end{bmatrix} \\begin{bmatrix} W \_q ^{11} & W \_q ^{12} \\\\ W \_q ^{21} & W \_q ^{22}\\end{bmatrix} \\begin{bmatrix} x \_m ^1 \\\\ x \_m ^2\\end{bmatrix} $$ Those familiar with linear algebra might recognize the [rotation matrix](https://en.wikipedia.org/wiki/Rotation%5Fmatrix?ref=lucasvandijk.nl). The input thus gets rotated, with the rotation angle being dependent on the token position $m$. We generalize to embedding dimensions greater than two by expanding matrix $\\textbf{R}$ as follows: $$ \\textbf{R}^d = \\begin{bmatrix} \\cos m\\theta\_1 & -\\sin m\\theta \_1 & 0 & 0 & \\ldots & 0 & 0\\\\ \\sin m\\theta \_1 & \\cos m\\theta \_1 & 0 & 0 & \\ldots & 0 & 0\\\\ 0 & 0 & \\cos m\\theta \_2 & -\\sin m\\theta \_2 & \\ldots & 0 & 0\\\\ 0 & 0 & \\sin m\\theta \_2 & \\cos m\\theta \_2 & \\ldots & 0 & 0\\\\ \\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots\\\\ 0 & 0 & 0 & 0 & \\ldots & \\cos m\\theta \_{d/2} & -\\sin m\\theta \_{d/2}\\\\ 0 & 0 & 0 & 0 & \\ldots & \\sin m\\theta \_{d/2} & \\cos m\\theta \_{d/2}\\end{bmatrix}$$ Pairs of features get rotated with varying base frequencies $\\theta \_i$. $\\theta\_i$ is defined similarly as in the absolute positional embeddings, with lower embedding dimensions rotating more quickly than higher embedding dimensions: $$ \\theta \_i = \\frac{1}{L^{2i/d}} $$ The benefit of using rotated vectors is that the angle between two position-adjusted vectors at $m$ and $n$ will only depend on their relative position, not on the absolute positions $m$ and $n$. ![Animation showing how the rotation angle between two vectors only depends on their relative position.](https://lucasvandijk.nl/content/images/2025/08/RoPEScene_ManimCE_v0.19.0.gif) The angle between two vectors depends on their relative position. ## Interpreting input features as complex numbers to derive the RoPE equation How was RoPE's rotation equation obtained? To explain its derivation, let's recall that the attention weights are computed by taking the softmax of the inner product between *query* and *key* vectors. In the case of *self-attention,* the query and key vectors are derived from different tokens in the sequence: $$ a\_{nm} = \\frac{\\exp{\\langle \\textbf{q}\_n^{T}, \\textbf{k}\_m\\rangle}}{\\Sigma\_{m'=1}^N \\exp{\\langle \\textbf{q}\_n^{T}, \\textbf{k}\_{m'}\\rangle}} $$ When using absolute positional encodings, the query and key vectors will be derived from the position-adjusted vector $\\textbf{x}\_n'$. For RoPE, we aim to find an alternative function $f$ to encode positional information, such that the inner product between two position-adjusted vectors only depends on the relative position between two tokens. Mathematically, the following should hold: $$ \\langle f(\\textbf{q}, m), f(\\textbf{k}, n)\\rangle = g(\\textbf{q}, \\textbf{k}, n - m) \\tag{1}$$ To derive such a function, we turn to the complex domain. A *d-*dimensional vector $\\textbf{v}\_n \\in \\mathbb{R}^d$ can be interpreted as a complex (*d/*2)-dimensional vector $\\textbf{v}\_n' \\in \\mathbb{C}^{d/2}$. For example: $$ \\textbf{v}\_n = \\begin{bmatrix} v\_1\\\\ v\_2\\\\ v\_3\\\\ v\_4\\end{bmatrix} \\rightarrow \\textbf{v}\_n' = \\begin{bmatrix}v\_1 + v\_2i\\\\ v\_3 + v\_4i\\end{bmatrix}$$ A complex number $a + bi$ can alternatively be written in polar form $R e^{i\\theta}$, where $R$ and $\\theta$ are called the *radial* and *angular* component of a complex number, respectively, and which are defined as $R = |a + bi| = \\sqrt{a^2 + b^2}$ and $\\theta = \\arctan(b/a)$. The standard inner product on $\\mathbb{C}^d$ is defined as $\\langle \\textbf{u}, \\textbf{v}\\rangle = \\overline{\\textbf{u}}^T\\textbf{v}$, i.e., taking the complex conjugate and transpose of $\\textbf{u}$ before multiplication with $\\textbf{v}$. Using this definition, and writing the complex numbers returned by $f$ and $g$ in polar form, we can obtain the following: $$ \\overline{R\_f(\\textbf{q}, m) e^{i\\Theta\_f(\\textbf{q}, m)}} \\cdot R\_f(\\textbf{k}, n) e^ {i\\Theta\_f(\\textbf{k}, n)} \\\\= R\_g(\\textbf{q}, \\textbf{k}, n-m) e^ {i\\Theta\_g(\\textbf{q}, \\textbf{k}, n-m)} \\tag{2}$$ Here we have explicitly written the inner product between two complex numbers on the left-hand side, equating it to a complex number on the right-hand side, representing $g(\\textbf{q}, \\textbf{k}, n-m)$. Additionally, we have defined functions $R\_f(\\textbf{q}, m)$, $\\Theta\_f(\\textbf{q}, m)$, $R\_g(\\textbf{q}, \\textbf{k}, n-m)$, and $\\Theta\_g(\\textbf{q}, \\textbf{k}, n-m)$, representing the radial and angular components of the complex numbers returned by $f$ and $g$, respectively. The conjugate of a complex number in polar form $Re^{i\\theta}$ is $Re ^{-i\\theta}$. With this definition, we can further simplify the left-hand side of Equation 2 by combining the product of two complex numbers: $$ R\_f(\\textbf{q}, m) R\_f(\\textbf{k}, n) e^{i(\\Theta\_f(\\textbf{k}, n) - \\Theta\_f(\\textbf{q}, m))} = R\_g(\\textbf{q}, \\textbf{k}, n-m) e^{i\\Theta\_g(\\textbf{k}, \\textbf{q}, n-m)} \\tag{3}$$ Equating the radial and angular components in Equation 3, we obtain: $$ R\_f(\\textbf{q}, m) R\_f(\\textbf{k}, n) = R\_g(\\textbf{q}, \\textbf{k}, n-m) \\\\ \\Theta\_f(\\textbf{k}, n) - \\Theta \_f(\\textbf{q}, m) = \\Theta\_g(\\textbf{q}, \\textbf{k}, n-m) $$ To make functions $R\_f$ and $R\_g$ more concrete, we start with a few simple initial conditions. Let's require that a token at the start of a sequence is unmodified by our positional embedding: $f(\\textbf{v}, 0) = \\textbf{v}$. Here $\\textbf{v}$ can refer to $\\textbf{q}$ or $\\textbf{k}$. Additionally, we define $\\textbf{v} = ||\\textbf{v}|| e^{i\\theta\_v} = R\_f(\\textbf{v}, 0) e^{i\\Theta\_f(\\textbf{v}, 0)}$, i.e., $R\_f(\\textbf{v}, 0) = ||\\textbf{v}||$, and $\\Theta\_f(\\textbf{v}, 0) = \\theta\_v$, where $\\theta\_v$ is the angular component of $\\textbf{v}$, i.e., the angular component of the reinterpreted values of $\\textbf{q}$ or $\\textbf{k}$ as a complex number. Now, let's assess the case when $n=m$: $$ R\_f(\\textbf{q}, m) R\_f(\\textbf{q}, m) = R\_g(\\textbf{q}, \\textbf{k}, 0) = R\_f(\\textbf{q}, 0) R\_f(\\textbf{k}, 0) = ||\\textbf{q}|| ||\\textbf{k}|| \\tag{4}$$ Here, we have applied Equation 3 to obtain the result $R\_g(\\textbf{q}, \\textbf{k}, 0)$, and applied it another time to obtain $R\_g(\\textbf{q}, \\textbf{k}, 0) = R \_f(\\textbf{q}, 0) R \_f(\\textbf{k}, 0)$. We use our initial condition definitions to obtain the final result. Note that Equation 4 should hold for all token positions $m$; the final result, however, doesn't depend on $m$. Thus the general solution for $R\_f(\\textbf{x}, m) = ||\\textbf{x}||$. What can we infer for the radial component when $n=m$? Applying Equation 3 similarly to what we did for the radial component, we obtain: $$ \\Theta\_f(\\textbf{k}, m) - \\Theta \_f(\\textbf{q}, m) = \\Theta \_g(\\textbf{q}, \\textbf{k}, 0) = \\Theta \_f(\\textbf{k}, 0) - \\Theta \_f(\\textbf{q}, 0) = \\theta\_k - \\theta\_q \\tag{5}$$ Rearranging terms of the first and last components of Equation 5, we obtain: $$ \\Theta\_f(\\textbf{k}, m) - \\theta\_k = \\Theta\_f(\\textbf{q}, m) - \\theta\_q $$ This holds for all $\\textbf{q}, \\textbf{k}$ and $m$, which implies that the change in the angular component of the position-adjusted complex number is independent of the $\\textbf{q}$ or $\\textbf{k}$. This insight gives us some freedom to define $\\Theta\_f(\\textbf{v}, m)$. For example, we can define it as follows: $$ \\Theta\_f(\\textbf{v}, m) = \\theta\_v + \\phi(m) $$ In other words, we simply add a position-dependent value $\\phi(m)$ to the original angular component $\\theta\_v$. But what would $\\phi(m)$ be? We can gain some insight by analyzing the difference between two subsequent positions, i.e., the case when $n = m +1$. $$ \\Theta\_f(\\textbf{k}, m+1) - \\Theta \_f(\\textbf{q}, m) = \\theta \_k + \\phi(m+1) - \\theta\_q - \\phi(m) \\\\ \\Theta \_f(\\textbf{k}, m+1) - \\Theta \_f(\\textbf{q}, m) + \\theta\_q - \\theta\_k = \\phi(m+1) - \\phi(m)\\\\ \\Theta \_g(\\textbf{q}, \\textbf{k}, 1) + \\theta\_q - \\theta\_k = \\phi(m+1) - \\phi(m)$$ Here we moved $\\theta\_q$ and $\\theta \_k$ to the left-hand side to obtain an expression for $\\phi(m+1) - \\phi(m)$. Furthermore, by applying Equation 3, we transformed the left-hand side to an expression *that does not depend on $m$.* This implies that the difference $\\phi(m+1) - \\phi(m)$ does not depend on $m$ either, and can be constant. We can define $\\phi(m)$ as a simple [arithmetic progression](https://en.wikipedia.org/wiki/Arithmetic%5Fprogression?ref=lucasvandijk.nl): $$ \\phi(m) = m\\theta + \\gamma$$ In this definition, $\\theta$ and $\\gamma$ are free to choose. For example, we can set $\\theta$ to an embedding-dimension-specific value, and set $\\gamma = 0$. Combining all pieces, we arrive at the following: $$ R\_f(\\textbf{v}, m) = ||\\textbf{v}||\\\\\\Theta \_f(\\textbf{v}, m) = \\theta \_v + m\\theta\\\\ f(\\textbf{v}, m) = ||\\textbf{v}||e^{i\\theta \_v + m\\theta} = \\textbf{v}e^{im\\theta}$$ Through our constraint on relative position and a few simple initial conditions, we were able to define a simple expression that encodes the sequence position into a query or key vector. Finally, using [Euler's formula](https://en.wikipedia.org/wiki/Euler%27s%5Fformula?ref=lucasvandijk.nl), we obtain the rotation matrix we discussed at the beginning of this post: $$ \\begin{aligned}e^{im\\theta} &= \\cos m\\theta + i\\sin m\\theta = \\cos m\\theta \\begin{bmatrix}1 & 0\\\\0 & 1\\end{bmatrix} + \\sin m\\theta \\begin{bmatrix} 0 & -1 \\\\1 & 0\\end{bmatrix} \\\\&= \\begin{bmatrix} \\cos m\\theta & -\\sin m\\theta\\\\ \\sin m\\theta & \\cos m\\theta\\end{bmatrix}\\end{aligned} $$ ## Conclusion Token order is critical for accurately interpreting natural language and biological sequences. By constraining the inner product between position-adjusted query and key vectors to be only dependent on their relative position, RoPE derives a simple transformation, rotating input vectors in the complex domain. This ensures the attention weights for two tokens are also influenced by their relative position. ## References Bishop, C. M., Bishop, H. & Cham, S. (ed.) (2023). *Deep Learning - Foundations and Concepts*. ISBN: 978-3-031-45468-4 Su, Jianlin, Yu Lu, Shengfeng Pan, Ahmed Murtadha, Bo Wen, and Yunfeng Liu. “RoFormer: Enhanced Transformer with Rotary Position Embedding.” arXiv:2104.09864\. Preprint, arXiv, November 8, 2023\. [https://doi.org/10.48550/arXiv.2104.09864](https://doi.org/10.48550/arXiv.2104.09864?ref=lucasvandijk.nl). Vaswani, Ashish, Noam Shazeer, Niki Parmar, et al. “Attention Is All You Need.” arXiv:1706.03762\. Preprint, arXiv, August 2, 2023\. [https://doi.org/10.48550/arXiv.1706.03762](https://doi.org/10.48550/arXiv.1706.03762?ref=lucasvandijk.nl). ### Fast and exact gap-affine partial order alignment with POASTA URL: https://lucasvandijk.nl/publications/fast-exact-gap-affine-partial-order-alignment-with-poasta/ Last updated: 2025-05-01T01:45:58.000Z Biology is full of sequences. An organism's hereditary information is encoded as a DNA molecule, and chains of amino acids fold into proteins, the cell's major functional molecules. Comparing such sequences is at the heart of many bioinformatic analyses. For example, it enables inferring the evolutionary history of a set of sequences (e.g., genes) or linking specific mutations to phenotypes. Such analyses require sequences to be *aligned,* i.e., we have matched the nucleotides or amino acids sharing an evolutionary ancestor among the different input sequences. To illustrate this, consider two seemingly different sequences `ATGCTTA` and `TGCAATCA`. We identify several shared nucleotides `(|)` by introducing *gaps* `(-)`and allowing *mismatches* `(*)` (Figure 1). ``` ATGC--TTA ||| |*| -TGCAATCA ``` ****Figure 1\.** Example alignment between two DNA sequences. Experts in their domain used to painstakingly create such alignments by hand. These days, we use algorithms to compute them. Pairwise alignment algorithms compute alignments between two sequences, while *multiple* sequence alignment (MSA) algorithms compute alignments for a set of sequences, e.g., for a set of genes. [Partial order alignment](https://academic.oup.com/bioinformatics/article/18/3/452/236691?ref=lucasvandijk.nl) (POA) is one type of MSA algorithm with applications in pangenome graph construction, *de novo* genome assembly, and many others. For example, the recently released [human pangenome reference](https://www.nature.com/articles/s41586-023-05896-x?ref=lucasvandijk.nl) was created with an algorithm that includes POA. POA can additionally be used to error-correct reads, e.g., as a preprocessing step before *de novo* genome assembly. POA represents the MSA as a *directed acyclic graph* (DAG) and works by iteratively computing sequence-to-graph alignments, updating the graph to include the new sequence after each alignment (Figure 2). The original algorithm introduced by Lee *et al.* is a dynamic programming approach, and the most recent implementation of this algorithm is called [SPOA](https://github.com/rvaser/spoa?ref=lucasvandijk.nl). SPOA utilizes SIMD instructions on modern CPUs, efficiently computing the dynamic programming matrix. ![](https://lucasvandijk.nl/content/images/2025/05/Screenshot-2025-04-30-at-9.37.05-PM.png) ****Figure 2\. (top)** A multiple sequence alignment represented as a directed acyclic graph. The alignment on the right represents the alignment of another sequence aligned to a path in the graph represented by CCGCAAAACGGG. ****(bottom)** The updated graph after incorporating the additional alignment in the top panel. POASTA is another approach to POA. It reformulates each sequence-to-graph alignment as a mathematically equivalent "shortest path" problem. Finding shortest paths is a common computer science problem, and many existing algorithms exist. Two classic algorithms are [Dijkstra's algorithm](https://en.wikipedia.org/wiki/Dijkstra%27s%5Falgorithm?ref=lucasvandijk.nl) and the [A\* algorithm](https://en.wikipedia.org/wiki/A%2A%5Fsearch%5Falgorithm?ref=lucasvandijk.nl). POASTA uses A\* as its main algorithmic framework and includes three algorithmic innovations to accelerate alignment: 1) a novel A\* heuristic that guides the algorithm more quickly to the target, 2) a depth-first extension component that exploits matches between the query sequence and the graph, and 3) it uses the graph topology to detect alignment states that will not be part of the optimal solution and skips those (Figure 3). These techniques substantially reduce the required computations to find the optimal alignment. ![](https://lucasvandijk.nl/content/images/2025/05/poasta-intro-a-2.png) ![](https://lucasvandijk.nl/content/images/2025/05/poasta-intro-b-2.png) ![](https://lucasvandijk.nl/content/images/2025/05/poasta-intro-c-3.png) ****Figure 3\.** POASTA accelerates alignment using three algorithmic innovations. ****(a)** POASTA introduces a "minimum gap cost" A\* heuristic. **(b)** POASTA exploits exact matches between the query sequence and the graph in a depth-first manner. ****(c)** POASTA uses "superbubbles" to detect and prune suboptimal alignment states. Besides accelerating alignment, it reduces the required amount of memory, enabling alignment of longer sequences than what was possible before. [In the paper](https://academic.oup.com/bioinformatics/article/41/1/btae757/7942505?ref=lucasvandijk.nl), we created MSAs of 342 \~1 megabase *Mycobacterium tuberculosis* sequences! POASTA is written in Rust and is available under the BSD-3 license. [GitHub - broadinstitute/poasta: Fast and exact gap-affine partial order alignmentFast and exact gap-affine partial order alignment. Contribute to broadinstitute/poasta development by creating an account on GitHub.![](https://lucasvandijk.nl/content/images/icon/pinned-octocat-093da3e6fa40-1.svg)GitHubbroadinstitute![](https://lucasvandijk.nl/content/images/thumbnail/poasta)](https://github.com/broadinstitute/poasta?ref=lucasvandijk.nl) ### Gut microbiome dysbiosis linked with recurrent UTIs URL: https://lucasvandijk.nl/publications/gut-microbiome-dysbiosis-linked-with-recurrent-utis/ Last updated: 2025-04-15T23:23:17.000Z Urinary tract infections (UTIs) are common, painful, and a large burden on healthcare worldwide. They frequently become recurrent, with some getting an infection multiple times a year. Furthermore, treatment becomes more challenging with increasing antibiotic resistance among the bacteria that cause UTIs. It is thus important to understand what drives UTI recurrence and devise treatment strategies that do not require antibiotics. *Escherichia coli* is the most common bacterial species causing UTIs, and the gut is a known reservoir of *uropathogenic* (UPEC) strains. Less understood, however, is the role of the entire gut microbiome in facilitating recurrence and whether any immunological differences play a role. [Longitudinal multi-omics analyses link gut microbiome dysbiosis with recurrent urinary tract infections in women - Nature MicrobiologyMulti-omics analyses of faecal, urine and blood samples from women with and without recurrent urinary tract infections reveal that gut dysbiosis and differential immune responses may play a role in risk of infection via the gut–bladder axis.![](https://lucasvandijk.nl/content/images/icon/apple-touch-icon-f39cb19454.png)NatureColin J. Worby![](https://lucasvandijk.nl/content/images/thumbnail/41564_2022_1107_Fig1_HTML.png)](https://www.nature.com/articles/s41564-022-01107-x?ref=lucasvandijk.nl) In this paper, we studied women with a history of recurrent UTIs for a year, collecting monthly blood, stool, and urine samples, and additional samples when diagnosed with a UTI. This enabled us to characterize the dynamics of the gut microbiome and explore its relationship with recurrent UTIs. Blood samples additionally provided a snapshot of their immunological state. We compared women in our UTI cohort with a matched cohort of healthy women. Our main findings include (1) the gut microbiome of women with UTIs are significantly less diverse, with fewer bacterial species present than typically seen in healthy women and showing more characteristics of low-level inflammation, (2) immunological biomarkers suggested a distinct immunological state, (3) *E. coli* strains frequently transmitted between the gut and the bladder in both the healthy cohort and the UTI cohort, though healthy women did not exhibit any UTI symptoms, and (4) the UTI-causing strain was rarely cleared from the gut after antibiotic treatment. The latter two findings were enabled by StrainGE, a software tool we specifically designed to identify and characterize low-abundance strains in complex microbial communities. It can untangle same-species strain mixtures and identify strain-specific genetic variants, enabling the tracking of strains across samples. [StrainGE: Strain Genome ExplorerStrainGE is a toolkit for tracking and characterizing low-abundance strains in complex microbial communities. It enables detailed insights into the bacterial strain-level diversity of whole metagenomic sequencing samples.![](https://lucasvandijk.nl/content/images/icon/favicon-2.ico)Lucas van DijkLucas van Dijk![](https://lucasvandijk.nl/content/images/thumbnail/strainge-cover-2.png)](https://lucasvandijk.nl/publications/strainge-strain-genome-explorer/) ### StrainGE: Strain Genome Explorer URL: https://lucasvandijk.nl/publications/strainge-strain-genome-explorer/ Last updated: 2025-04-16T18:04:42.000Z We humans each carry thousands of bacterial species. They inhabit our skin, nose, gut, and several other sites. Many species protect us from pathogens, help digest food, or train our immune system. Other species, however, can cause life-threatening infections. Even within the same bacterial species, enormous differences in pathogenicity exist. For example, almost every human harbors the bacteria Escherichia coli in their gut microbiome without any issues. Some strains, however, can cause severe diarrhea or urinary tract infections. Because of this phenotypic diversity, it is important to know what specific strain is present when analyzing microbial communities. Improved insights into the strain-level diversity of complex microbial communities will strengthen our understanding of their role in human health. For example, by comparing the genomes of strains, we could identify genetic factors distinguishing pathogenic and non-pathogenic strains. [StrainGE: a toolkit to track and characterize low-abundance strains in complex microbial communities - Genome BiologyHuman-associated microbial communities comprise not only complex mixtures of bacterial species, but also mixtures of conspecific strains, the implications of which are mostly unknown since strain level dynamics are underexplored due to the difficulties of studying them. We introduce the Strain Genome Explorer (StrainGE) toolkit, which deconvolves strain mixtures and characterizes component strains at the nucleotide level from short-read metagenomic sequencing with higher sensitivity and resolution than other tools. StrainGE is able to identify strains at 0.1x coverage and detect variants for multiple conspecific strains within a sample from coverages as low as 0.5x.![](https://lucasvandijk.nl/content/images/icon/apple-touch-icon-582ef1d0f5.png)BioMed CentralLucas R. van Dijk![](https://lucasvandijk.nl/content/images/thumbnail/13059_2022_2630_Fig1_HTML.png)](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02630-0?ref=lucasvandijk.nl) This is where StrainGE comes in. StrainGE is specifically designed to identify and characterize low-abundance strains using metagenomic sequencing data from a community. Metagenomic sequencing data represents sequenced DNA fragments from all community members, which is a challenge for characterizing low-abundance strains. For example, *E. coli* typically represents only 1% of a healthy human gut microbiome; thus, only a tiny fraction of sequenced reads will originate from *E. coli* strains. StrainGE overcomes this challenge by comparing the read data to a database of reference genomes, such as genomes available in public databases. It first detects which reads in the sample likely originate from the species of interest and then compares the reads to the references in the database, reporting those that look the most similar to the strain(s) in the sample. If multiple strains of the same species are present, it will report multiple reference genomes. StrainGE's reported references serve as a basis for further, more detailed characterization of the sample strains. Since the reported references are unlikely to be the same as the strains in the sample, StrainGE identifies strain-specific genetic variants by mapping sample reads to the references. StrainGE analyzes the read alignment pileups to search for evidence of different alleles compared to the reference. ![](https://lucasvandijk.nl/content/images/2025/04/strainge-overview.webp) Overview of the StrainGE algorithm. StrainGE was instrumental in characterizing the *E. coli* strain-level dynamics in the gut microbiomes of women with recurrent urinary tract infections (UTIs). We found that the UTI-causing strain was rarely cleared from the gut after antibiotics and found unexpected similarities with a healthy control group. More about this study can be found at the link below. [Gut microbiome dysbiosis linked with recurrent UTIsMore than half of the women in the US get a urinary tract infection (UTI) in their lifetime, which frequently becomes recurrent. In this paper, we investigated the role of the gut microbiome in facilitating recurrence.![](https://lucasvandijk.nl/content/images/icon/favicon-7.ico)Lucas van DijkLucas van Dijk![](https://lucasvandijk.nl/content/images/thumbnail/umb-nmicro-2-1.jpg)](https://lucasvandijk.nl/publications/gut-microbiome-dysbiosis-linked-with-recurrent-utis/) ### Using locality sensitive hashing to compactly represent k-mers URL: https://lucasvandijk.nl/2018/03/using-locality-sensitive-hashing-to-compactly-represent-k-mers/ Last updated: 2025-04-14T21:33:14.000Z Comparing (DNA) sequences is one of the core tasks in bioinformatics, and the classic approach is to align these sequences. This is, however, a relatively slow process and not always computationally feasible, especially if you want to compare more than two DNA sequences. An alternative approach is to compare sequences based on their *k-mer profiles*. A *k-mer* of a string $S$ is defined as any substring of $S$ of length $k$. For example, the DNA sequence `AGCGTATCGATTCA` has the following k-mers if $k=6$: ``` AGCGTATCGATTCA -------------- AGCGTA GCGTAT CGTATC ... GATTCA ``` As you can see, obtaining all *k-mers* is easy: slide a window of size $k$ along your sequence, yielding a k-mer at each position. A sequence of length $L$ has $L−k+1$ k-mers. A common task is to count how often each k-mer occurs and compare genomes based on these counts. The main idea is that [similar genomes have similar *k-mer* counts](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1319-7?ref=lucasvandijk.nl). When dealing with the scale of genomes, storing counts for all these different *k-mers* can take up quite a lot of memory. First, the number of distinct *k-mers* grows exponentially with the length of $k$. In the case of DNA sequences, our alphabet size is 4: A, C, T, G. Therefore, there are $4^k$ possibilities of length $k$. The value of $k$ depends on your application and organism, but values ranging from 5 to 32 are common. Next, think how we would store the *k-mer* itself. We could store each letter as ASCII character, requiring 8-bits per character. However, this is a bit wasteful because we only have four characters in DNA. An optimisation would be to use 2 bits per character: A=00, C=01, T=10, G=11\. This would allow us to store a *k-mer* of length 32 in a 64-bit integer. Still, this may not be good enough. I’ve seen cases for $k=23$ where it went up to more than 100 GB, and that’s quite a lot of memory even if you have access to a decent compute cluster. This post will explain a technique described in the paper by [Cleary *et al*. for to reduce the memory consumption for storing *k-mers*](https://www.nature.com/articles/nbt.3329?ref=lucasvandijk.nl). The main insight is that we often don’t need the *exact* count of each *k-mer*, and can take some shortcuts by missing some *k-mers*. Because of the exponential number of different *k-mers* and because genomes are often large, missing a few *k-mers* will not have a huge impact. Furthermore, when dealing with whole genome sequencing datasets, we must also deal with sequencing errors and expect some *k-mers* to be false. In a lot of cases, using approximate *k-mer* counts is appropriate. We start by transforming each *k-mer* to a complex vector, using the following encoding: $$A=1,T=−1,C=i,G=i$$ A k-mer is now a vector of length $k$, where each element represents a DNA base as a complex number. A *k-mer* is now a point in a k-dimensional space. Notice that similar *k-mers* will be “neighbours” in this k-dimensional space. This is not an efficient encoding, but we can apply a few smart operations to convert this to a smaller integer. The idea is to divide our k-dimensional space into a lot of bins. This can be done as follows: generate a random complex vector of length k, which can then be interpreted as the normal vector of a plane through our k-dimensional space. This plane divides our space in half: certain k-mers lie on the “left” side, while other k-mers lie on the “right” side of this plane. See the figure below for an example. If we note a zero when the *k-mer* lies on the left and a one if on the right, we obtain a bit value. ![](https://lucasvandijk.nl/content/images/2025/04/kmer-space-split.png) ****Figure 1\.** A few (poorly drawn) examples of splitting a k-dimensional space in half by a hyperplane. Repeat this n times, and you obtain n bit values. For example, we could store the result from drawing 32 random hyperplanes in a 32-bit integer. Note that we split the space in half each time we draw a random hyperplane, creating $2^n$ bins in our $k$-dimensional space. It is possible, however, that two distinct *k-mers* end up in the same bin and thus have the same 32-bit integer value (if $n$=32). We can actually calculate the probability of such an event. Recall that our *k-mers* are just ordinary complex vectors and that similar *k-mers* are neighbours in the k-dimensional space. We can compute the cosine of the angle between the two vectors as follows: $$ \\phi = \\cos(\\theta) = \\frac{\\textbf{u}\\cdot\\textbf{v}}{||\\textbf{u}||||\\textbf{v}||} $$ Here **u** and **v** are two *k-mers* in our complex k-dimensional space. The probability of these two *k-mers* being in the same bin, is equal to the probability that *no* random hyperplane cuts through the angle of **u** and **v**. This is visualised in the figure below. ![](https://lucasvandijk.nl/content/images/2025/04/kmer-neighbours.png) ****Figure 2\.** The probability that a hyperplane cuts through the angle of two similar k-mers is lower than with two very dissimilar k-mers. We end up with the following equation: $$ P(h(\\textbf{u}) = h(\\textbf{v})) = 1 - \\frac{\\arccos(\\phi)}{\\pi} $$ Here, $h(u)$ represents the *hash* value of a *k-mer* **u**, or the n-bit integer by generating $n$ random hyperplanes. This is the reason why this could be used for *approximate* counting of *k-mers*, because some *k-mers* may be mistaken for another. By tuning $k$ and $n$ you can try to minimize the impact on your application. The effect is even further minimized because similar *k-mers* have a higher probability of ending up in the same bin than dissimilar *k-mers*, while reducing the number of bits you need to store the *k-mer*. ### Visualizing the height of the Netherlands URL: https://lucasvandijk.nl/2018/02/visualizing-the-height-of-the-netherlands/ Last updated: 2025-04-15T20:45:12.000Z On this date 65 years ago, February 1st 1953, the Netherlands experienced its greatest flood till date, the [North Sea Flood of 1953](https://en.wikipedia.org/wiki/North%5FSea%5Fflood%5Fof%5F1953?ref=lucasvandijk.nl). This is still one of the biggest disasters the Netherlands has ever experienced, with thousands of casualties and lots of people who lost their homes. The Netherlands earns its name because large parts of the country lie below sea level. To make sure our country doesn’t flood, we have built lots of barriers, dams and dykes to keep the water out, and we want to prevent anything like the flood of 1953 from happening ever again. In the beginning of this year, several of these barriers and dams were put to the test when a heavy storm reached the Netherlands which resulted in very high water levels. Our five biggest dams and barriers needed to be closed at the same time, a first since their construction. We can ask ourselves the question whether this will happen more often now that [sea levels are rising due to global warming](https://climate.nasa.gov/vital-signs/sea-level/?ref=lucasvandijk.nl). A higher base line sea level increases the chance for even higher water levels when it storms. To get an idea what areas would be affected the most by a possible flood, we have created a visualisation project that shows the height of the Netherlands in comparison to the sea level. ![](https://lucasvandijk.nl/content/images/2025/04/maeslantkering.jpg) The Maeslantkering in closed state. Photo: [Rijkswaterstaat](https://beeldbank.rws.nl/?ref=lucasvandijk.nl). ## Things to check out - The Y-axis shows deviation from “[Amsterdam Ordnance Datum](https://en.wikipedia.org/wiki/Amsterdam%5FOrdnance%5FDatum?ref=lucasvandijk.nl)“ - You can zoom on both the map and the elevation profile (by selecting a region in the bottom profile). - If you hover your mouse over the elevation profile, a red dot will show you the location on the map - The highest point in our country is the [Vaalserberg](https://nl.wikipedia.org/wiki/Vaalserberg?ref=lucasvandijk.nl), with a whopping 322 m! It’s in the most southern part of the Netherlands. Some people, however, may claim it’s not really part of the Netherlands because of their ugly accents ;) - If you zoom in on our shorelines, you can see our dykes ## Elevation profiles of the Netherlands ## Technical details ### Quick links - Data source: [http://ahn.nl](http://ahn.nl/?ref=lucasvandijk.nl) - We use AHN version 2 and the dataset with a resolution of 5 meter - Written in Python and CoffeeScript, and available on Github: [https://github.com/lrvdijk/nl-height](https://github.com/lrvdijk/nl-height?ref=lucasvandijk.nl) ### Algorithm overview 1. The publicly available height data is split across multiple “raster data” files. The AHN website provides a “ahn\_units.shp” shapefile, that shows which part of the Netherlands is covered by which *raster data* file, and thus shows which file to download. 2. We import this ahn\_units shapefile in a PostgreSQL database with the PostGIS extension enabled. A script downloads and extracts all available raster data files. 3. When all data is downloaded, we import the raster data in our PostgreSQL database. 4. To generate an elevation profile, we draw a line horizontally across the Netherlands and query which cells in our raster data hit our line. Each cell contains the height at its midpoint. Our query is based on this blogpost: [http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html](http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html?ref=lucasvandijk.nl) 5. When there’s water, no data exists in our database. We apply some post-processing such that our visualisation does not get weird artefacts. 6. We store the data in a JSON file, and use the d3.js library for the visualisation. ### Drawing arbitrary shapes with OpenGL points URL: https://lucasvandijk.nl/2015/06/drawing-arbitrary-shapes-with-opengl-points/ Last updated: 2025-04-15T20:45:03.000Z Part of my Google Summer of Code project involves porting several arrow heads from [Glumpy](https://github.com/glumpy/glumpy?ref=lucasvandijk.nl) to [Vispy](https://vispy.org/?ref=lucasvandijk.nl). I also want to make a slight change to them: the arrow heads in Glumpy include an arrow body, I want to remove that to make sure you can put an arrow head on every type of line you want. Making a change like that requires that you understand how those shapes are drawn. And for someone without a background in computer graphics this took some thorough investigation of the code and the techniques used. This article is aimed at people like me: good enough programming skills and linear algebra knowledge, but almost no former experience with OpenGL or computer graphics in general. ## Implicit surfaces The basic principle behind drawing these 2D shapes is called implicit surfaces. It relies on an arbitrary shape function that returns the distance to your shape surface or boundary from a given point in your image. This is visualized in Fig. 1. ![](https://lucasvandijk.nl/content/images/2025/04/fig1-implicit-surfaces.png) ****Figure 1\.** Distances to the boundary of a circle. The distances any shape function should return are highlighted in red. To actually be able to draw a shape, we need to distinguish whether a point lies inside the shape or not. We make the arbitrary decision that a negative distance lies inside a shape, and a positive distance means that the point lies outside the shape. ## Distance functions for a few basic shapes The distance functions defined below have one requirement: the center point of the shape has the coordinate (0, 0). ### Circle For a circle these distances are easily calculated: $$ d(\\textbf{x}) = ||\\textbf{x}|| − r $$ Where: - $\\textbf{x}$: Vector to the point in consideration. - $r$: The radius of the circle. If the point lies within the circle, the length of the vector towards that point is smaller than the radius, making the distance automatically negative. ### Square Consider Fig. 2. ![](https://lucasvandijk.nl/content/images/2025/04/fig2-distance-square.png) ****Figure 2\.** Visualizing the distances to the border of a square. A square is a nice symmetric figure, so we can take the absolute values of the coordinates. Then, we only need to consider the smaller highlighted square (light blue). The distance to the boundary of the square is then: $$ d(\\textbf{x}) = \\text{max}(|x\_1|, |x\_2|) - \\frac{s}{2}$$ Where: - $\\textbf{x}$: Vector towards the point in consideration - $|x\_1|,|x\_2|$ are the absolute values of the first and second element of the vector (the *x* and *y* coordinates). - $s$: The size of the square. Using the max function we sort of select to which boundary the distance will be calculated. Then we substract the size of the smaller square (highlighted with light blue). The resulting distance is then negative if the point lies within the square, and positive otherwise. ## Combining shapes Combining multiple simple shapes is often useful to make more complex shapes. To do this, we introduce some basic operations on the returned distances of a simple shape. ### Inverse We arbitrarily decided to say that the distance is negative when a point lies within the shape. To get the inverse of a shape, we simply need to negate each distance value. $$\\forall x, y: \\neg S(x, y) = -S(x, y)$$ ### Union The union of two shapes can be retrieved by using the min function on both distance functions. $$ \\forall x, y : U(x, y) = \\min(S\_1(x, y), S\_2(x, y)) $$ Remember that the distance value is negative when the point lies in the shape. The lowest value will win here, so if one of those distances is negative (the point belongs to at least one shape), it will return the negative value. Thus combining both shapes to a single one. ### Difference The difference of two shapes contains all points in $S\_1$ excluding the points in $S\_2$. This is defined as follows: $$\\forall x, y : D(x, y) = \\max(S\_1(x, y), -S\_2(x, y))$$ Consider the example where $S\_1(x,y)=−2$ and $S\_2(x,y)=−1$. In short, the current point (x,y) belongs to both S1 and S2\. Using the above function for the difference, the value from $S\_2$ gets negated: $−S\_2(x,y)=1$. Due to the max function, this value will win (it’s higher than -2), and therefore, it will not be part of the new shape. This is precisely what we want because we want all points that are part of $S\_1$ but not part of $S\_2$. ### Intersection The intersection of two shapes contains all points that are both part of $S\_1$ and $S\_2$. It is defined as follows: $$\\forall x, y: I(x, y) = \\max(S\_1(x, y), S\_2(x, y))$$ A point will be part of the new shape if and only if both distances are negative. If one distance is positive, the max function will return this value, and a positive value means it is not part of the shape. This results in a shape that includes only points that are both part of $S\_1$ and $S\_2$. ## OpenGL implementation So, how do we translate these principles into working code? Let's first introduce some basic OpenGL concepts before we present the shader code. ### Shaders and drawing modes I will not go too deep in the basics of OpenGL, but a modern OpenGL pipeline consists of multiple *shaders*, small programs you can write yourself. At the very minimum you’ll need a *vertex shader* and a *fragment shader*, which determine where the main “drawing points” will be positioned and the color of the individual pixels respectively. The pipeline is illustrated in Fig. 3. ![](https://lucasvandijk.nl/content/images/2025/04/fig3-gl-pipeline.png) ****Figure 3\.** OpenGL pipeline illustrated. Courtesy of [Nicolas Rougier.](http://www.labri.fr/perso/nrougier/teaching/opengl/?ref=lucasvandijk.nl) You can define your own attributes for each vertex, for example, the position, color, or orientation. OpenGL has several modes for generating the actual primitives. These are illustrated in Fig. 4. ![](https://lucasvandijk.nl/content/images/2025/04/fig4-gl-primitives.png) ****Figure 4\.** OpenGL primitive generation modes. Courtesy of [Nicolas Rougier.](http://www.labri.fr/perso/nrougier/teaching/opengl/?ref=lucasvandijk.nl) For a more in-depth OpenGL introduction, I recommend [this tutorial](http://www.labri.fr/perso/nrougier/teaching/opengl/?ref=lucasvandijk.nl), [Anton’s OpenGL tutorials](http://antongerdelan.net/opengl/?ref=lucasvandijk.nl), or [opengl-tutorial.org](http://opengl-tutorial.org/?ref=lucasvandijk.nl). ### General 2D shape shaders For the 2D shapes we want to draw, we’ll use the points drawing mode. OpenGL allows you to specify the point size, and the fragment shader will be called for each pixel in the point. Let’s check the vertex shader where we position our vertices and configure the point size. #### Vertex shader Listing 1: The vertex shader code for our 2D shapes ```glsl // Uniforms // ------------------------------------ uniform float antialias; uniform mat4 ortho; // Attributes // ------------------------------------ attribute vec2 position; attribute float size; attribute vec4 color; attribute float rotation; attribute float linewidth; // Varyings // ------------------------------------ varying float v_size; varying vec4 v_color; varying vec2 v_rotation; varying float v_antialias; varying float v_linewidth; // Main // ------------------------------------ void main (void) { v_size = size; v_linewidth = linewidth; v_antialias = antialias; v_color = color; v_rotation = vec2(cos(rotation), sin(rotation)); gl_Position = ortho * vec4(position, 0, 1); gl_PointSize = M_SQRT2 * size + 2.0 * (linewidth + 1.5*antialias); } ``` We first define some *uniforms*, *attributes*, and *varyings*. Uniforms are variables which are the same for each vertex. Attributes are variables defined per vertex, and with varyings we can pass some data to the next steps in the pipeline (for example, the fragment shader). Each vertex has a position where our 2D shape will be drawn. The matrix `ortho` is used for the proper projection of the vertex to your screen. We will not explain this in-depth, but if you want to know more please refer to [this tutorial](http://www.opengl-tutorial.org/beginners-tutorials/tutorial-3-matrices/?ref=lucasvandijk.nl) on opengl-tutorial.org. Our shapes are larger than one pixel, so we need to change `gl_PointSize`. Our shapes also have a size attached to them, but for the point size we scale this size with $\\sqrt{2}$ (ignore the extra size for antialias en linewidth for now). We do this because our shapes can be rotated. To fit a rotated square of size *x* in another square, we need a bigger square of size $x\\sqrt{2}$ (I hope you remember Pythagoras). This is visualized in Fig. 5. ![](https://lucasvandijk.nl/content/images/2025/04/fig5-point-rotation.png) ****Figure 5.** Rotation of a square. Furthermore, we pass along some variables to the next steps in the pipeline (size, line width, antialias, color, rotation). Note we create a direction vector for the rotation from the given rotation in radians. #### Fragment shader Listing 2: The fragment shader code ```glsl // Varyings // ------------------------------------ varying float v_size; varying vec4 v_color; varying vec2 v_rotation; varying float v_antialias; varying float v_linewidth; // Main // ------------------------------------ void main() { vec2 P = gl_PointCoord.xy - vec2(0.5,0.5); P = vec2(v_rotation.x*P.x - v_rotation.y*P.y, v_rotation.y*P.x + v_rotation.x*P.y) * v_size; float size = v_size/M_SQRT2; float distance = shape_func(P, size); gl_FragColor = filled(distance, v_linewidth, v_antialias, v_color); } ``` Note that we have the same varying definitions as in the vertex shader. These contain values as passed from the vertex shader. Also note the usage of the built-in variable `gl_PointCoord`. We specified in the vertex shader the size of our point in pixels, and for each pixel in this point the fragment shader gets called. The `gl_PointCoord` contains the coordinates **inside the point**. The *x* and *y* attributes from `gl_PointCoord` range from 0.0 to 1.0, where (0, 0) is the bottom left corner of the point, and (1, 1) is the top right corner of the point. There are several operations applied to these coordinates: **Step 1\.** First we substract $\\begin{bmatrix}0.5 \\\\ 0.5 \\end{bmatrix}$. This ensures the origin is in the center of the point because the distance functions we defined earlier in this article require that. **Step 2.** Next**,** we apply a rotation transformation to the point. Remember the transformation matrix is: $$\\begin{bmatrix}\\cos(\\theta) & -\\sin(\\theta)\\\\ \\sin(\\theta) & \\cos(\\theta) \\end{bmatrix} $$ If you look closely at the code you see that the vector $\\textbf{P}$ gets multiplied with this matrix. $$\\begin{bmatrix}nx \\\\ ny\\end{bmatrix} = \\begin{bmatrix}x \\\\ y \\end{bmatrix} \\begin{bmatrix} \\cos(\\theta) & -\\sin(\\theta)\\\\ \\sin(\\theta) & \\cos(\\theta) \\end{bmatrix} = \\begin{bmatrix} x \\cos(\\theta) - y \\sin(\\theta) \\\\ x\\sin(\\theta) + y \\cos(\\theta) \\end{bmatrix} $$ **Step 3.** We also scale the coordinates with `v_size`. These transformations are visualized in Fig. 6. ![](https://lucasvandijk.nl/content/images/2025/04/fig6-transformations.png) ****Figure 6\.** Transformation of the coordinates within a point. The green region in the bottom of fig. 6 is our canvas for drawing our shape. This is done by `shape_func()`, any function that returns the distance to the boundary of a shape as explained earlier in this article. The `filled()` function determines the color for the current pixel determined by the returned distance of `shape_func()`. Simply put: if the returned distance is negative, it returns a color (because it’s part of the shape). Otherwise, it makes the current pixel transparent. It also applies some anti-aliasing techniques which I don’t know the details of, so we will not cover this in-depth. ## Example: curved arrows To conclude this article, we will examine the distance function of a curved arrowhead. A curved arrowhead can be constructed using the inverse of the union of three circles. This is visualized in Fig. 7, and the accompanying shader code can be seen in lst. 3. ![](https://lucasvandijk.nl/content/images/2025/04/fig7-curved-arrows.png) ****Figure 7\.** Construction of a curved arrow head. Listing 3: GLSL function to the distance of an arrow ```glsl /** * Computes the signed distance to an curved arrow * * Parameters: * ----------- * texcoord : Point to compute distance to * size : size of the arrow head in pixels * * Return: * ------- * Signed distance to the arrow * */ float arrow_curved(vec2 texcoord, float size) { vec2 start = -vec2(size/2, 0.0); vec2 end = +vec2(size/2, 0.0); float height = 0.5; vec2 p1 = start + size*vec2(0, -height); vec2 p2 = start + size*vec2(0, +height); vec2 p3 = end; // Head : 3 circles vec2 c1 = circle_from_2_points(p1, p3, 6.0*size).zw; float d1 = length(texcoord - c1) - 6*size; vec2 c2 = circle_from_2_points(p2, p3, 6.0*size).xy; float d2 = length(texcoord - c2) - 6*size; vec2 c3 = circle_from_2_points(p1, p2, 3.0*size).xy; float d3 = length(texcoord - c3) - 3*size; return -min(d3, min(d1,d2)); } ``` We first define the arrow corner points as `p1`, `p2` and `p3`. We use a helper function which calculates the center point of a circle through two points with a given radius `r`. The distance to these circles are then easily calculated, and we use the min function to get the union of these three circles. Our arrow head is exactly the area *not* covered by these circles (see fig. 7), and therefore we return the inverse of this union. ### Controlling a servo with an AVR microcontroller URL: https://lucasvandijk.nl/2012/10/controlling-a-servo-with-an-avr-microcontroller/ Last updated: 2025-08-12T21:53:27.000Z Servos are small motors that can be precisely controlled and often power the movements of small robotic arms. There are many tutorials on how to control them with an Arduino, but fewer tutorials use only a bare AVR chip. In this tutorial, we’ll be using the ATTiny44, a small and cheap microprocessor that also contains a 16-bit timer, which will make our lives a bit easier, as I'll explain later. Servos can be precisely rotated a specific number of degrees, depending on the pulse width you feed them with the microcontroller. They can also be used as motors to drive a wheel, though you’ll need special ‘continuous rotation’ servos. You’ll often find them in RC cars. So, let's get started and see how you control a servo! ## Pulse width modulation To move the servo, you need to send it a pulse. Depending on the duration of the pulse, the servo positions itself in a fixed position. This position is unique to the pulse duration. In other words, no matter what the current position of the servo is, it always rotates to the same position for a certain pulse duration. The following pulse durations result in the following positions: - 1ms: 90 degrees to the left - 1.5ms: Center position - 2ms: 90 degrees to the right For most servos, the maximum frequency is around 50 Hz, which means a time period of 20 ms. Of course, the servo is not limited to just 90-degree angles; the angle is proportional to the pulse duration, with a minimum of 0.7 ms and a maximum of 2.3 ms. ## Implementing it with an AVR microcontroller Let's start with the circuit, which is simple. ![Circuit with an AVR ATTiny44 to control a servo.](https://lucasvandijk.nl/content/images/2025/04/servoschema.jpg) Circuit with an AVR ATTiny44 to control a servo In this case, I use a 9V battery to power my AVR and 4 AA batteries (total of 6V) to power my servo. It’s better to use a different power source for your servos because they can draw high currents when they’re rotating, which could possibly trigger the AVR reset from the resulting voltage drop. I use the LM7805 to transform the 9V to 5V for the AVR, some capacitors to remove any AC components of the voltage, and a 4.7k resistor as a pull-up for the reset pin. ### AVR source code ```c /* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*- */ /* * main.c - Auto-generated by Anjuta's Makefile project wizard * */ #ifndef F_CPU #define F_CPU 1000000L #endif #include #include #include #define SERVO_PWM_TOP 10000UL // 50 Hz PWM #define SERVO_LEFT 500UL // Capture value to position the servo to the left (1.0ms) #define SERVO_CENTER 750UL // Capture value to position the servo in the center (pulse width 1.5ms) #define SERVO_RIGHT 1000UL // Capture value to turn the servo (pulse width 2ms) int main(void) { // Setup ports and timers DDRA = 0xFF; // All output PORTA = 0; // Configure timer/counter1 as phase and frequency PWM mode TCNT1 = 0; TCCR1A = (1 << COM1A1); TCCR1B = (1 << WGM13) | (1 << CS10); ICR1 = SERVO_PWM_TOP; OCR1A = SERVO_LEFT; while(1) { } } ``` So, what are we doing here? We will use the AVR's built-in PWM hardware to generate the desired pulse. The ATTiny44 contains a 16-bit timer, which has a special mode called ‘Frequency and phase correct PWM’. Grab the ATTiny44 datasheet and read about it. To summarize what it does, you provide a TOP value and a compare value. The AVR will then count to the TOP value, and after it has reached that, it counts back to zero. When the counter has the same value as the compare value in OCR1A, something will happen: in the case of up counting, it will reset the pin value of OCA1 to 0, and when down counting, it will set the pin value to 1\. This behaviour can be changed by setting the right bits in the TCCR1A/B/C registers. Read the datasheet for more information. To visualize a bit of what’s happening, here’s a drawing: ![Value of the timer, and the match events.](https://lucasvandijk.nl/content/images/2025/04/servo-timer-match.jpg) Value of the timer, and the match events. Some things you can conclude from the above drawing and description: The value in ICR1 defines the frequency, the value in OCR1A defines the duty cycle. ### Calculating the values for ICR1 and OCR1A So why does the above code has 10000 as value for ICR1, and 500/750/1000 for OCR1A? With a few simple calculations, you can get the values you need: 1. Check at which clockfrequency your AVR runs, in our case it’s 1 MHz. 2. Because it runs on 1 MHz, if we had an unlimited number of bits for our counter, it would count to 1000000 in one second. 3. Almost all servos operate at a frequency of 50 Hz, which means a time period of 20 ms. 4. You know this: 1000000 for one second, so 1000000 \* (20\*10^-3) = 20000 counts for 20 ms. 5. In phase and frequency correct mode, it counts up AND down, so we divide the above number by 2, which results in 10000, the value for ICR1\. This is why a 16-bit timer is useful: 10000 easily fits in 16 bits and definitely not in 8 bits. The same can be done to position the servo at the left (a pulse of 1 ms), or at the right (a pulse of 2 ms). **Please note that the above code will constantly pulse every 20 ms.** If the servo is already in the right position, it will not move. And it’s probably better to disable the timer, when you’ve sent the pulse. That’s left as an exercise for the reader. ;) ## Troubleshooting - In general, it’s better to use a separate power source for the servos to avoid triggering the AVR's reset. But remember: the grounds of each power source should be connected! - If you want your servo to rotate continuously, remember to buy one capable of doing that. Most servos can’t rotate a full 360 degrees. However, there are mods to modify a non-continuous rotation servo to one that can rotate the full 360 degrees. You can find a lot of tutorials on Google. And it’s pretty easy. Thanks for reading, and I hope you enjoyed the article! ### Multithreading with C++11: protecting data URL: https://lucasvandijk.nl/2012/06/multithreading-with-c-11-protecting-data/ Last updated: 2025-04-12T20:43:59.000Z Welcome to this second part in a series of articles about multithreading with C++11\. In the previous part, I briefly explained what a thread is and how to create one with the new C++ thread library. This time, we will be writing a lot more code, so open up your favourite IDE if you want to try the examples while you’re reading. In the previous article, we also saw that sometimes the output wasn’t completely right when running multiple threads simultaneously. Today, we’ll see some other problems with sharing a resource between threads and, of course, provide some solutions to these problems. *This is an article part of a series about multithreading with C++11, the other parts are listed below:* - Part 1: [Introduction to threads with C++11](https://lucasvandijk.nl/introduction-to-threads-with-c-11/) - Part 2: Protecting your data with multiple threads ## Lost data Let’s start with a simple example. Consider the following program: ```cpp #include #include #include #include using std::thread; using std::vector; using std::cout; using std::endl; class Incrementer { private: int counter; public: Incrementer() : counter{0} { }; void operator()() { for(int i = 0; i < 100000; i++) { this->counter++; } } int getCounter() const { return this->counter; } }; int main() { // Create the threads which will each do some counting vector threads; Incrementer counter; threads.push_back(thread(std::ref(counter))); threads.push_back(thread(std::ref(counter))); threads.push_back(thread(std::ref(counter))); for(auto &t : threads) { t.join(); } cout << counter.getCounter() << endl; return 0; } ``` The purpose of this program is to count to 300000\. Some smartass programmer wanted to optimize the counting and created three threads, each adding 100000 times one to a shared variable `counter`. Let’s walk a bit through the code. We create a new class called `Incrementer` which holds a private variable `counter`. The constructor is straightforward, simply initializing `counter` by setting it to zero. What follows is an operator overloading function, and in this case `operator()`. This means that each object of this class can be called as a simple function. Usually, you would call a method on an object like this: `object.fooMethod()`, but now, you can actually *call* the object, like this: `object()`. This is convenient because now we can pass the whole object to the thread class and, within the operator overloading function, use all the advantages of a class. The last method of the class is `getCounter`: a simple getter for the `counter` variable. Then, we have the `main()` function; here, we similarly create threads as described in the previous article. A few differences, though: we now create an object of the class `Incrementer`, and we pass it to the threads. Note that we use `std::ref` here, to pass a reference of the object, instead of passing a copy to the thread. So, let’s see what this program produces and if this smartass programmer is brilliant or stupid. Compile the program using GCC 4.7 or higher, or Clang 3.1 or higher. In the case of GCC, use the following command: `g++ -std=c++11 -lpthread -o threading_example main.cpp` And voilà, the output: ```bash [lucas@lucas-desktop src]$ ./threading_example 218141 [lucas@lucas-desktop src]$ ./threading_example 208079 [lucas@lucas-desktop src]$ ./threading_example 100000 [lucas@lucas-desktop src]$ ./threading_example 202426 [lucas@lucas-desktop src]$ ./threading_example 172209 ``` But wait a second! That smartass programmer wasn’t so smart after all! The program doesn’t count to 300,000; one run only reached 100,000! Why is this happening? Well, as simple as ‘increment by one’ sounds, it requires multiple instructions for the processor to increment a variable: ```asm movl counter(%rip), %eax addl $1, %eax movl %eax, counter(%rip) ``` What happens is this: 1. We have the current counter value stored in memory, and we load that value into the EAX register 2. We add one to the value in the EAX register 3. We store the value in EAX back in the original memory location I hear you think, “Ok, nice information, but how does this solve my counting problem, smartass?” Remember from the previous article that threads shared the processor if there’s only one core. So at some point, one thread is happy because its instructions are executed for a while, but then the big boss operating system says, “Ok, time’s up, back in the line!” Then, another thread will be executed for a while. And when it’s the turn of the original thread again, he starts executing where he was left. So, can you guess what happens when the operating system decides to switch to another thread while the original thread is in the middle of incrementing a variable's value by one? Well, for example, this: | Thread 1 | Thread 2 | Counter | Explanation | | --------------- | --------------- | ------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | | %eax <- counter | nothing | 1 | The current value of the counter is now stored in EAX | | %eax + 1 | nothing | 1 | The value in EAX is now 2, but the counter hasn't been updated yet | | nothing | %eax <- counter | 1 | Hey, a switch to the other thread, and because the counter hasn't been updated, it loads the old value | | %eax -> counter | nothing | 2 | And we're back at thread 1\. When a thread switches, the state of the registers (for example EAX), are saved, and when we switch back, they're restored in their original state. Because EAX was two the last time thread one got interrupted, that value will be written in the counter variable. | | nothing | %eax + 1 | 2 | And we're back at thread 2\. The last time thread got interrupted, EAX had the value 1\. | | nothing | %eax -> counter | 2 | And because EAX was one (and is now two because of the previous instruction), the counter remains 2. | ## Hey, it’s occupied! The solution is to make sure only one thread can access the shared variable at the same time. This can be done using the [std::mutex](http://en.cppreference.com/w/cpp/thread/mutex?ref=lucasvandijk.nl) class. Visualize it as a sort of toilet: when you go inside, you lock it, do your stuff, and then you unlock it. Anyone who wants to use the bathroom has to wait before you’re ready. A convenient feature of a mutex is that the operating system makes sure the locking and unlocking operations are indivisible. This means that the thread will not be interrupted when it’s trying to lock or unlock a mutex. When a thread locks or unlocks a mutex, this operation will be finished before the operating system switches threads. And the best thing is, when you try to lock a mutex, but some other thread has already locked it, you’ll have to wait. But the operating system keeps track of which threads wait on which mutex. The blocked thread will go into a “blocked on m” state, meaning the operating system won’t give that thread any processor time until the mutex becomes unlocked. This means no wasted CPU cycles. If multiple threads are waiting, it depends on the operating system which thread will be the ‘winner.’ General purpose operating systems like Linux and Windows use a First In First Out system; on realtime operating systems, it’s priority-based. Let’s modify the above code to make sure counting works as expected. ```cpp #include #include #include #include #include using std::thread; using std::vector; using std::cout; using std::endl; using std::mutex; class Incrementer { private: int counter; mutex m; public: Incrementer() : counter{0} { }; void operator()() { for(int i = 0; i < 100000; i++) { this->m.lock(); this->counter++; this->m.unlock(); } } int getCounter() const { return this->counter; } }; int main() { // Create the threads which will each do some counting vector threads; Incrementer counter; threads.push_back(thread(std::ref(counter))); threads.push_back(thread(std::ref(counter))); threads.push_back(thread(std::ref(counter))); for(auto &t : threads) { t.join(); } cout << counter.getCounter() << endl; return 0; } ``` Note the changes: we included the `mutex` header file and added a member `m` to our class, with type [mutex](http://en.cppreference.com/w/cpp/thread/mutex?ref=lucasvandijk.nl), the standard mutex class in C++11\. In the method `operator()()`, we lock the mutex just before incrementing the counter, and afterwards, we unlock it again. When we run the program, the output is correct: ```bash [lucas@lucas-desktop src]$ ./threading_example 300000 [lucas@lucas-desktop src]$ ./threading_example 300000 ``` And as always in computer science, there’s no free lunch. Using mutexes will considerably slow down your program, but that’s better than an incorrect program. ## Heimdall, guard us against exceptions When incrementing a variable by one, the chance of an exception being raised is not particularly high, but with more complex code, it’s definitely possible. The above code is not really exception-safe. When an exception occurs, the mutex will still be locked while the function is already finished. To make sure the mutex is unlocked when an exception is thrown, we could use the following code: ```cpp for(int i = 0; i < 100000; i++) { this->m.lock(); try { this->counter++; this->m.unlock(); } catch(...) { this->m.unlock(); throw; } } ``` But that’s an awful lot of code for just locking and unlocking a mutex. Luckily, there’s a nice, simple solution for that: the [std::lock\_guard](http://en.cppreference.com/w/cpp/thread/lock%5Fguard?ref=lucasvandijk.nl) class. The `lock_guard` class is really simple: it locks the given mutex on creation, and unlocks the mutex when the lock is destroyed (for example, at the end of a function scope). Modifying the above code again results in: ```cpp void operator()() { for(int i = 0; i < 100000; i++) { lock_guard lock(this->m); // The lock has been created now, and immediatly locks the mutex this->counter++; // This is the end of the for-loop scope, and the lock will be // destroyed, and in the destructor of the lock, it will // unlock the mutex } } ``` This code is also exception-safe because when an exception occurs, the lock's destructor will still be called, resulting in an unlocked mutex. Remember, you can create your temporary scopes in the following way: ```cpp void long_function() { // some long code // Just a pair of curly braces { // Temp scope, create lock lock_guard lock(this->m); // do some stuff // Close the scope, so the guard will unlock the mutex } } ``` ## Closing I promised in the previous article to explain why the printing of “Hello World” and “Parallel World” sometimes went wrong. I hope you understand now that sharing resources between threads without synchronization can cause problems. I will be a bit lazy and redirect you guys [to this stackoverflow question](http://stackoverflow.com/questions/6374264/is-cout-synchronized-thread-safe?ref=lucasvandijk.nl). To fix the output of cout, create another mutex, and each thread should lock that mutex before sending something to cout and, of course, unlock it when it’s done. I hope you enjoyed the article. Next time, we’ll cover condition variables, another widely used technique to synchronize threads. ### Introduction to threads with C++11 URL: https://lucasvandijk.nl/2012/05/introduction-to-threads-with-c-11/ Last updated: 2025-04-06T17:34:02.000Z [The free lunch is over](http://www.gotw.ca/publications/concurrency-ddj.htm?ref=lucasvandijk.nl). The time when our complex algorithm ran extremely slow on a computer but ran extremely fast a few years later because the processor speeds exploded is gone. The trend with current processors is to add more cores rather than increase the clock frequency. As a programmer, you should be aware of this. Of course, processors will always perform a bit better each year, but the growth in performance is slowing down. Many programs can benefit the most by using multiple threads because of today’s multicore processors. In this article, I’ll briefly explain what a thread is and how you can create them with the new threading library in C++11\. I’m planning to write multiple articles about this topic, each going more in-depth. *This is part of a series of articles about multithreading with C++11, the other parts are listed below:* - Part 1: Introduction to threads with C++11 - Part 2: [Protecting your data when using multiple threads](https://lucasvandijk.nl/2012/06/multithreading-with-c-11-protecting-data/) ## What is a thread? On most general-purpose computers, we run many processes alongside each other. But let us ask the question, what is a process? Well, a simple definition is a running program. But what is a program, then? Another simple definition is a list of instructions the processor needs to execute. Assume we have one single-core processor. How can all these programs run beside each other? After all, a single-core processor can perform only one instruction at a time. Well, that’s handled by the operating system. To share the processor with all these processes, it gives one process a little bit of time to perform some instructions, then goes to another process which gets a little bit of time, and so on. Some processes have higher priority than others and will get more processor time. This is called scheduling and is a subject on itself, and I won’t cover it much more in this article. We can have threads within processes. It’s based on the same principle described above: on a single-core processor, each thread is given some processor time. If we have a non-threading program, we start at the `main()` function, and the program is finished when we reach the end of the `main()` function. Fun fact: When you run this program, the operating system actually creates a new thread to run the `main()` function in, which we call the ‘main thread’ (creative name, isn’t it?). You can create new threads from the main thread, which should run simultaneously with the main thread. This means that besides the ‘codepath’ `main()` till the end of `main()` we now also have another codepath, from the entry point of the new thread till the end of the thread. The thread's entry point is often the start of a function other than `main()`, and the end of the thread means the end of that function. But remember, on a single-core processor, it’s not really simultaneous; it shares the processor with the other threads. We discussed single-core processors in the above text, but what about today’s multicore processors? Well, this means we can do multiple things at the same time. If your processor has two cores, then two threads can run simultaneously. If your processor has N cores, then your processor can run N threads at the same time. And this is why today multithreading is so important: we can actually do things in parallel. ### Difference between processes and threads Threads and processes both have the same purpose: running specific tasks simultaneously. Why do they both exist? Well, there are some differences, summarized below. #### Properties of a process - The stack of a process is safe. - Each process has its own memory, which other programs can’t alter. (There are probably ways to do it, but in normal circumstances, it can’t be done.) - Because each process has its own memory, the memory is **safe**, but inter-process communication is *slow* - Switching from one process to another is heavy: processor-cache flush, memory management unit TLB flush. #### Properties of a thread - The stack of a thread is safe - Each thread shares the same memory within the same process - Shared memory is **unsafe**, but inter-process communication is *fast* - Switching from one thread to another is *fast* (no flushes) ## Defining threads with C++11 Fine, fine, let’s start coding. C++11's new threading library is really nice, and makes creating a new thread easy. Consider the example below; we’ll walk through the code afterward. ```cpp #include #include #include #include #include using std::string; using std::thread; using std::vector; using std::cout; using std::endl; /** * Prints the given string `num_times` times, with `time_between` miliseconds * between each print. */ void printer(string text, int num_times, int time_between) { for(int i = 0; i < num_times; i++) { cout << text << endl; // Wait a moment before next print std::chrono::milliseconds duration(time_between); std::this_thread::sleep_for(duration); } } int main() { // Create the threads which will each print something vector printers; printers.push_back(thread(printer, "Hello", 30, 10)); printers.push_back(thread(printer, "World", 40, 15)); printers.push_back(thread(printer, "Parallel World", 40, 20)); for(auto &p : printers) { p.join(); } return 0; } ``` Of course, we first include some library files. The `chrono` and `thread` libraries are the most important here. Both are new in C++11\. The `chrono` library provides some nice timing capabilities, and the `thread` library the classes for creating threads. Then, we define a new function called `printer`. This is a straightforward function that just prints the given text a number of times, with a given time in between. You can see that the `chrono` library provides a nice and clean way to define the time between each print. The `main` function is a bit more interesting; here, we create the other threads. C++11 provides a class `thread`. As you can see, we create three objects of this class, and put them all in a vector. The `thread` constructor accepts a function as the first argument, and the rest of the arguments will be passed to the given function. Because the `vector` implements the iterator API, we can use the new C++11 range based for syntax, and on each `thread` object we call the method `join`. This ensures the calling thread will now wait for the joined thread to exit before it finishes itself (in this case, the main thread, and thus the whole program). In our case, each printer thread must exit before the main thread finishes. ### Running the program When we compile and run this program, we get the following output: ``` HelloWorld Parallel World Hello World Parallel World Hello World Hello Parallel World Hello World Hello Parallel World World Hello Hello World Parallel World Hello World Hello Parallel World World Parallel World World World Parallel World Parallel World Parallel World Parallel World Parallel World Parallel World Parallel World Parallel World ``` So yeah, we can see that our printer threads run beside each other, and not sequential. There’s one thing to note here, although we have the following code: `cout << text << endl;` In the first line of the above output, the ‘Hello’ and ‘World’ are not on seperate lines. Something is not entirely right, but the why and how to fix it will be covered in my next article.