In a select few corn fields after a harvest, the unwanted stalks, leaves, and cobs are fed into a storage container on the back of a long-haul truck. These agricultural wastes won’t be hitting the road in their loose, fluffy state though. The storage container doesn’t start empty. It contains a complex chemical machine fine-tuned to efficiently convert the corn detritus into three chemical compounds: biochar, bio-oil, and syngas. The biochar, a solid, is scattered back onto the farm fields to improve field fertility. The syngas is, in part, used as fuel to keep the high temperature reaction going. The bio-oil is destined to be injected deep underground. Bizarrely, people will pay for it to be put there, not to save for later use, but with the promise that it will stay there forever.

This whole operation happens under the purview of Charm Industrial, a leading participant in a growing market of players with a shared goal: suck carbon dioxide directly out of the air and store it away for good. This is carbon dioxide removal, or CDR for short.

The machinery acting as this conceptual vacuum cleaner for the sky is about as diverse as you can imagine for an industrial process. Charm’s bio-oil production is just one of many. Another method starts with trays of powdered rock and ends with cement production. A third involves growing kelp in the open ocean and encouraging it to eventually sink to the bottom.

The market is small for now, but ambitious. To meet the scale of needs for carbon dioxide removal outlined in the Intergovernmental Panel on Climate Change’s (IPCC) 2021 report, our rate of global removal will need to double roughly every 21 months for the next 27 years ¹. We’ll need a market with hundreds of billions of dollars in transactions per year. We’ll need an entirely new global industry on the scale of steel or cement.

This is the rate we need assuming we succeed in dramatically reducing emissions and dramatically scaling up natural climate solutions.

So we’ve got some challenges ahead.

Since leaving my job as a software engineer at Figma at the end of 2021, I’ve been learning about frontier problems within climate, and seeing where software can play a role. Unfortunately, my exploration was perilously broad. Climate tech encapsulates everything from predicting the cost of connecting new solar to the grid, to chemically transforming refrigerants, to improved cow feed to reduce methane from their burps, to novel steel manufacturing strategies. Each entails a totally different market, with different stakeholders, and different problems.

Without choosing a specific problem domain, I found myself forever skimming across the surface: learning small details here and there, but never learning enough to evaluate real hypotheses about where I, personally, could apply leverage.

So lately I’ve been focusing on a narrow but important slice: carbon dioxide removal. I’d like to share what I’ve learned about the CDR market through reading and a couple dozen conversations with founders, investors, engineers, and scientists. This article is for anyone curious about CDR, but has a slight tilt towards areas where software is playing a role.

Here’s a preview of the path we’ll take:

Before we get started, let’s locate ourselves in the space, and select our lens.

Orienting ourselves in the broader context

Talking about carbon dioxide removal is a bit of a mess.

I’m afraid readers might leave with misguided perspectives like “new technology is the solution to climate change”, or worse “carbon dioxide removal is the solution to climate change”. New technology broadly, and carbon dioxide removal specifically, are squarely in the camp of “necessary, but not sufficient”.

We’ll be exploring a sub-branch of a sub-branch of what’s needed to stabilize our climate. But remember, we need the whole conceptual tree to win this fight, and watchful arborists to tend to each branch and ensure its growth.

As for real trees, they will not be our focus today. Nor will filters placed atop smoke stacks. We’ll be focused on permanent removal.

Three crucial jobs

If you examine this graph describing an optimistic path to net-zero emissions, the three colors (, , ) indicate three crucial jobs needed to reach a stable climate over the coming decades:

1. Reduce new emissions (shrink the grey ). This is where the vast majority of the work lies in fighting climate change. This job requires transitioning our electrical grid to renewable energy, electrifying our transportation and heating, and decarbonizing industrial processes like cement and steel production. This is also where technology like filters on smoke stacks, or filters placed on the exhaust system of semi-trucks belong. Some of these solutions, like replacing a coal power plant with a solar farm plus storage, might be permanent in the sense that the emissions are permanently reduced. But solutions in this camp can’t remove emissions from the atmosphere that happened long ago. So things in this camp are maybe permanent, but not removal.
2. Use nature to remove carbon from the atmosphere (grow the green ). This is where high-quality tree planting initiatives play a role, for example TIST and Terraformation. These do remove carbon dioxide directly from the atmosphere, even if it was emitted long ago. But we can’t be confident how long carbon held in trees will stay there. Forest fires, lumber production, and disease can all result in the carbon returning to the atmosphere. It is removal, but it’s not confidently permanent.
3. Artificially remove carbon from the atmosphere (grow the blue ). There are a few reasons we need to do this on top of the other two jobs. First, we don’t have credible plans for avoiding certain kinds of emissions, like the partial evaporation of fertilizer from agricultural fields². That’s why you never see the grey section go away completely. Second, we want to return our atmosphere to pre-industrial levels of CO2 concentration, not just stop concentrations from rising. That requires removing more carbon-dioxide than we’re emitting, which is why you see the black line representing net emissions eventually going into the negative.

Today we’ll only be discussing the third job. I didn’t choose to focus on CDR because it’s the only important job of these three, or even the most important of these three. We need all three. But you can’t develop real insight into any of them if you don’t focus for a while.

The lens of the day

There are many lenses to use when trying to understand CO2 removal. We could ask how we developed consensus for how much we need by 2050, or what spurred the initial research identifying strategies for CO2 removal, or compare different methods by their energy and land area requirements. Those all fascinating angles! But they’re not the focus of this article. Here, our goal is to build intuition for what the industry looks like in 2023.

In terms of big buckets, we’ll be focusing on three:

1. Suppliers (): the companies doing the actual removal, and the partners they need to make it happen
2. Capital (): where suppliers get the money to do their job
3. Trust (): how the suppliers and capital build and maintain consensus about how much carbon is being removed, and who paid for it

We’ll build up the following map, but do so gradually, so don’t worry about imbibing its contents in a single gulp. We’ll use concrete examples of organizations playing these various roles, and aim to be illustrative, rather than exhaustive.

With boundaries drawn around the subject of the day, we can start to examine what’s inside those boundaries. Let’s start with the suppliers.

Suppliers: the companies doing the main thing

To be a supplier of permanent CDR, you have two crucial jobs:

1. Capture: how do you pull CO2 from the atmosphere?
2. Sequestration: where do you put the captured CO2, and how do you ensure it stays there for 1000+ years?

Let’s see how some leading suppliers answer these two questions:

The problems these companies are facing in scaling removal are varied, but I heard a few problems repeated. One was a problem I never needed to think about for software: permitting.

Permitting: for anything that touches air, water, or dirt

To deploy at scale, removal suppliers need to build industrial infrastructure. Whenever you want to build a big physical anything, permitting law comes into play. Discovering which permits you need, applying for them, and engaging in back-and-forth with their regulating body can be a massive time sink for suppliers.

When I visited Charm Industrial, CEO Peter Reinhardt explained that they don’t deploy industrial infrastructure in California, despite company HQ being in San Francisco. Getting permits in California can take years. By contrast, in some more politically right-leaning states, it takes weeks. They’d love to deploy in California, but the timelines make it prohibitive.

When I asked Heirloom’s CEO Shashank Samala about the problems that keep him up at night, permitting was on his list:

Different jurisdictions have varying rules and requirements for permits such as grading, electrical, building, construction, environmental and other permits. It would be immensely beneficial to have a comprehensive and up-to-date database containing all the necessary information for each permit, in every jurisdiction including contact information and where exceptions can be made. Oftentimes, the challenge lies in the fact that permitting teams are understaffed. We had issues with staff turnover / lost knowledge midway through given how long these permitting timelines are. It is astonishing how outdated some of these procedures can be. For instance, in California, physical signatures on laminated drawings are still required for civil drawing approvals.

A few early-stage startups have begun working in this space. Blumen Systems and Paces are both working on automating discovery of permits and siting intelligence for climate tech companies.

The permitting problem is painful enough that CDR suppliers may design their deployment strategy around it. For some companies, this means deploying in states or countries with less arduous permitting processes. Another strategy is to woo the right industrial partners.

Industrial partners: piggy-back on a massive industry

Many CDR providers form symbiotic relationships with industry interests. The industrial partners can accelerate CDR scaling by allowing the CDR supplier to operate on their land and use their permits. The CDR companies can save their partners money by turning their partners’ waste products into tools for CDR.

For example:

• Charm Industrial has agricultural and forestry partners to turn waste streams (corn stover and mechanical forest thinning respectively) into feedstock for its pyrolyzers. The agricultural partners get biochar, a soil-health enhancing material, distributed onto their fields in the process.
• Ebb Carbon partners with desalination plants, helping them turn the waste hyper-saline solution extracted from sea water into a chemical input for Ebb. By installing on-site, Ebb gains access to a system for returning large volumes of water to the ocean.
• Lithos has mining and agricultural partners. They turn basalt (a waste product of the mining industry) into soil health enhancing minerals for their agricultural partners. The agricultural partners provide them the large land area they need for spreading the basalt to maximize its reactivity.

With lab-tested CDR methodologies and industrial partners in hand, suppliers can build out a plan to scale. To put that plan into action, they need cash, and lots of it.

Capital: show me the money

The CDR suppliers have bills to pay, things to buy, and employees that need to eat. CDR suppliers’ capital comes predominately from selling carbon removal credits, venture capital, and research grants.

Carbon removal credits

Suppliers want to sell the removal of CO2 as a service. To make this a purchasable quantity, they parcel up tons of CO2 removed into units sometimes called “carbon removal credits”. The credit will show, at a minimum, how much CO2 was removed, who did the removal, and when the removal happened. Buyers of these credits can use this information as an evidence of fulfilling public facing climate commitments.

So who are these buyers? Currently, nobody has a regulatory requirement to buy carbon removal ³, so anyone doing so is doing it voluntarily.

The industry is targeting $100/ton of CO2 removed, and 3.8 billion tons/year by 2050, which would make for a $380 billion/year market. At that scale, private-sector voluntary markets are going to be woefully inadequate, but we need to bootstrap a market to get there. This requires some buyers to act in economically irrational ways.

The corporate buyers

Surprisingly, there are a number of companies willing to do it anyway. Some are driven by a sense of moral obligation to build a better future. Some may expect this market to happen with or without them, so they’d like a seat at the table to shape future legislation. Whatever the motive, it’s been enough to kickstart a sizeable market. According to cdr.fyi, over 3 million tons in sales have been made as of May 2023, for a total transaction volume of $365 million USD since the industry’s inception.

One of the earliest buyers was Stripe ⁴, who announced in 2019 that they’d spend at least $1 million USD/year on carbon removal. In 2020, they fulfilled that promise by making purchases from Climeworks, Project Vesta, Charm Industrial, and Carbon Cure ⁵. For 3 of these companies, Stripe was their first customer.

Nan Ransohoff, Stripe’s current Head of Climate, now also leads Frontier, a coalition of Stripe, Alphabet, Shopify, Meta, and McKinsey which have collectively committed $1 billion towards purchasing carbon removal between 2022 and 2030 ⁶. This kind of advanced market commitment is crucial for building financial security for suppliers, which in turn provides security for investors knowing that these companies will have revenue to bolster their growth.

On May 23, 2023, JPMorgan Chase signed contracts to purchase $200 million in carbon removal via Climeworks, Charm Industrial, and a contribution to Frontier.

Microsoft explicitly committed to carbon negativity, which requires the purchase of carbon removal. They have a live dashboard of the removals they’ve purchased contracts for, including which vendors provided it.

The curators

Early in this process, all purchases were made by bilateral agreements between climate-conscious companies and carbon removal providers, with no market curators in the middle. This meant that buyers had to do due diligence from scratch, and providers had to spend a ton of time on sales. This is one of the inefficiencies that companies like Watershed and Patch are addressing.

The equivalent for individuals instead of corporations also exists: Wren and Commons provide subscriptions to individuals who want to offset their carbon footprint, and they buy carbon offsets and removals on your behalf. If you’re looking for a tax-deductible option, check out Terraset.

Pre-sale of credits for fledgling companies

The money available to carbon removal suppliers once they can deliver removal is encouraging! But the market, as it stands today, is heavily supply constrained. While 3 million tons have been purchased, only ~2% of those tons have actually been delivered. Many suppliers are sold out of commitments for the next several years. So at the moment, we need ways of accelerating our suppliers, and need ways of creating more suppliers. AirMiners is trying to accelerate the creation of CDR companies by helping them secure early funding through pre-sale of credits, discounted heavily based on the risk of credits never being delivered. Frontier also offers pre-purchase for fledgling CDR suppliers.

The more traditional route for pre-revenue companies to get funding is venture capital.

Venture capital

Venture capitalist funds raise money from wealthy individuals and funds-of-funds to invest in portfolios of high-risk companies. They assume that the majority of their portfolio companies fail, but that a few of them yield outsized returns. Making many bets in pure software is relatively cheap, because the cost to evaluate product-market fit is usually the salary of few software engineers for a year or two. CDR companies, by contrast, typically need millions of dollars just to buy equipment, even for lab-scale prototypes. Because CDR is fundamentally about scaling physical processes, both the cost and speed of scaling is much worse than software, making it less appealing to most venture capitalists.

There are a few firms, however, with enough conviction in the space to invest in multiple providers: ⁷

• Lowercarbon Capital invested in Charm Industrial, Heirloom, Running Tide, Verdox, Undo, and Noya. Lowercarbon has a $350 million fund dedicated to carbon removal.
• MCJ Collective invested in Charm Industrial, Heirloom, and Noya.
• Prime Impact Fund invested in Charm Industrial, Verdox, and Project Vesta.
• Breakthrough Energy Ventures invested in Heirloom, and Verdox.

Here’s what Yin Lu from MCJ Collective (a purely climate fund) said about how they think about investments into companies with physical assets as part of their portfolio:

Roughly half of our investments historically have been in pure software companies and we expect that trend to continue.

Roughly half of our investments have some sort of real-world component to them. […] Most of our physically oriented companies have created some form of new process such as carbon removal, chemical synthesis, food production, or the like. What you need to believe for this category is that the industries they are in will be undergoing fundamental change over the next decade due to the move toward decarbonization. This is a once in a lifetime transition that will see new companies emerge to capture large shares of existing markets and achieve scale that was otherwise not possible in previous innovation cycles.

Each time you raise venture capital, you sell part of the company. Founders looking for additional funding but wary of losing ownership will search for other sources of capital. The most common for early climate tech companies is grants.

Grants

For novel technology aligned with societal need, there are billions of dollars available in research grants.The government offers a variety of grants that are applicable to some subset of CDR companies:

• The Office of Clean Air Demonstrations has $3.5 billion available in grants to establish direct air capture hubs.
• The Advanced Research Projects Agency-Energy (ARPA-E) announced $45 million in funding for companies working on validation of ocean-based carbon dioxide removal and also provided funding for research into Direct Air Capture (DAC) technology, including a grant to Heirloom for $476,811.

Private organizations also offer grants:

• The Musk Foundation partnered with XPRIZE to offer a $100 million milestone-based prize for carbon dioxide removal providers.
• Additional Ventures created a $10 million research award to further research into ocean alkalinity enhancement.

Unlike venture capital, the money for grants is not exchanged for equity in the company receiving the grant. The granting agency gets nothing financial in return. This makes grants, in a sense, “free money”. But the application process can be time consuming, and grants commonly constrain how the money can be used. Even after grants are issued, they create additional work reporting progress back to the granting agency.

It’s common for CDR companies to have full-time on-staff grant writers or to hire consulting firms to help them fill grants to manage the paperwork load. Some companies skip grants altogether because of their time intensity. To help ease this process, companies like Streamline and Pioneer are leveraging modern AI (e.g. GPT-4) to gather requirements and fill in grant applications.

Even after grant approval, the time for the money to hit the company’s bank account can sometimes be upwards of a year. This is one place where loans can help.

Loans

If you’re quite confident you’ll have money in the future, but you need it now in order to help your company grow faster, you want some kind of loan. Enduring Planet builds financial infrastructure to handle two cases of this situation, specifically for climate entrepreneurs.

The first is when you’ve demonstrated a steady stream of revenue. This is unsurprisingly called “Revenue Based Financing”. The second, a “Climate Grant Advance”, provides money to climate companies who’ve already been approved for a grant, but haven’t received the money yet.

As the market matures and the revenue streams for CDR companies becomes less risky, traditional finance organizations will play a larger role here.

We’ve outlined two crucial blocks of the market: suppliers to remove carbon, and sources of money to make it happen. This system only works, however, if the folks with the money trust that their money has the intended effect.

Building trust: where, exactly, did my money go?

If you owned a grocery store and needed to source some tomatoes, you’d try a few suppliers, and evaluate them based on the price, quality, and reliable delivery. If the tomatoes taste like feet or always show up late, you’d switch suppliers. In the CDR market, customers are buying a promise that an invisible gas is moved from one place to another, where the destination is far out of sight from the customer. So not only is it hard for the customer to “taste” the quality of the product, it’s hard to even know when it’s arrived!

To make the system work, we need a group of third parties where their primary function is to facilitate trust between suppliers and sources of capital.

This starts with building a shared understanding of the core mechanism of removal.

Fundamental research: does the process remove carbon at all?

The discovery and evaluation of methods for carbon dioxide removal originate in fundamental research within universities and government labs. Through peer review and challenges from subsequent papers, we develop increasing confidence in the fundamental principle behind each method.

For example, the paper CO2 extraction from seawater using bipolar membrane electrodialysis (2012) was peer-reviewed and published in the journal of the Royal Society of Chemistry in their Energy & Environmental Science journal. The first author of the paper, Matthew Eisaman, went on to found Ebb Carbon to implement and scale the practice described in the paper.

Ambient weathering of magnesium oxide for CO2 removal from air (2020) was peer-reviewed and published in Nature Communications. The first author, Noah McQueen, is now co-founder of Heirloom, which is scaling this methodology.

This kind of research can serve as the catalyst to convince venture capitalists that a company concept isn’t totally nuts. As part of their investment memo for Heirloom, MCJ Collective cited this underlying science as a source of confidence:

[T]he basis of Heirloom’s technology stems from research led, in part, by Dr. Jennifer Wilcox who now serves at the U.S. Department of Energy. While there inevitably will be a need to bridge the gap between research performed at the lab bench and commercializing the technology in market, we are confident this risk is managed by the encouraging findings of the research and Heirloom’s iterative prototyping.

– “Our Investment in Heirloom”, MCJ Collective

These papers describe fundamental mechanisms and estimation protocols, but typically describe isolated steps at small scale. To understand climate impact, we need to look at the process holistically.

LCA: is it carbon negative?

Having a mechanism which removes carbon dioxide from the atmosphere is useless if you emit more carbon than you remove in the process. Emissions can come from e.g. energy used in manufacturing industrial equipment, operation of that equipment, and transportation. Determining the net removal is the domain of Life Cycle Analysis (LCA).

Absent relevant regulation, the process for establishing consensus on an LCA for a given company’s removal methodology is still up in the air, but typically involves peer-review from industry interests rather than academics. These LCAs are made publically available for scrutiny.

The framework Charm Industrial uses to evaluate its net emissions was written by Carbon Direct and EcoEngineers, and then reviewed by Charm Industrial, Lowercarbon Capital, Frontier, and others.

The framework Running Tide uses was written by Running Tide, then reviewed by Lowercarbon Capital, Stripe, Patch, and others.

Once the fundamental removal mechanism is confirmed, the carbon negativity of the entire process is quantified via LCA, and the supplier can run the process outside the lab, the last crucial step is providing an auditing chain for each action in the removal. This is the domain of Measurement, Reporting, and Verification (MRV).

MRV: what did my money, specifically, do?

To see what MRV looks like in practice, it’s instructive to look at Charm Industrial’s registry. On September 24, 2021, Block purchased 55 tons of CO2 removal from Charm Industrial via Watershed:

To present this information to buyers, you need:

• Sensors deployed in the field (in this case, probably a scale for bio-oil)
• A data collection pipeline to record this data reliably and accurately
• Software to convert those raw sensor values into emission and removal numbers using the appropriate method’s LCA process
• A way to build out a many-to-many relationship between deliveries and purchases and durably record those purchases to avoid double-counting

Most steps in Charm’s process lend themsleves well to direct measurement. If we trust Charm is telling the truth about the scale reading leading up to the bio-oil being injected deep underground, then we can trust the amount of removal being performed. This is also the case for direct air capture systems like Heirloom, where the pure CO2 stream captured can be directly measured.

There are still some sources of uncertainty in these processes, however. Using its CDR Verification Framework, CarbonPlan is working on mapping out the sources of uncertainty in calculating net emissions, and assigning “Verification Confidence Levels (VCLs)” ranging from 1 to 5 for different pathways.

CarbonPlan assigns “Biomass Carbon Removal and Storage”, the pathway used by Charm, a VCL of 3-5. To examine what the remaining sources of uncertainty are, check out the CarbonPlan website.

For CDR methods using the ocean or soil to sequester carbon, measurement is more difficult, which yields higher uncertainty. For example, for ocean-based methods to accurately estimate their net removal, they need:

• Access to third party meteorological and oceanographic data like weather predictions, ocean currents, seabed topography, and satellite imagery
• A much richer variety of sensors deployed into the ocean on buoys: humidity, barometric pressure, temperature, etc.
• More resilient systems for sensor data collection, since the sensors will need to be left unattended for long spans of time
• Computationally intensive ocean simulations to estimate the timeline and quantity of the removal

Unlike Charm or Heirloom’s processes, it’s infeasible to directly measure the amount of CO2 removed from the atmosphere. The removal takes place over a massive surface area of water, as the ocean itself absorbs CO2 from the atmosphere. Even if we all agree on all the sensor inputs, there’s a lot to debate about the methodology for turning those sensor inputs into the number of tons of CO2 removed.

These are the kinds of challenges that Running Tide and Ebb Carbon face. CarbonPlan assigns Running Tide’s original method, “Ocean Biomass Sinking (No Harvest)”, a VCL of 1-2, and Ebb Carbon’s method, “Ocean Alkalinity Enhancement (Electrochemical)”, a VCL of 3. We can raise these confidence levels over time by developing better measurement strategies through research and experimentation. This is presumably why, as noted in the grants section above, ARPA-E announced USD$45 million in funding to validate marine CDR methods.

If we can control uncertainty, ocean and soil-based systems hold a lot of promise. By piggy-backing off of natural processes already happening at massive scale, they tend to be easier to scale quickly with lower energy requirements and require less novel machinery ⁸.

At the moment, this MRV process is done in-house by the CDR suppliers. This both creates duplication of effort (e.g. every supplier building their own web interface to expose to buyers), and a perverse incentive to over-estimate removal. For now, the market is small enough (and difficult enough to enter as a supplier) that trusting companies to act in good faith is reasonable. As the market scales, however, this will become decreasingly true. So we’ll eventually need MRV-specific companies.

MRV service companies

A tiny minority of funding (especially venture capital) has gone to companies focused on MRV for CDR. But there are a few players emerging.

Isometric is building a multi-pathway registry and verification service. To minimize perverse incentives to fill their registry with as many credits as possible, they charge a per-ton verification fee to the buyers of carbon removal rather than accept per ton registration fees from suppliers. The registry portion of the platform would be similar to Charm Industrial’s registry, but supporting many companies and CDR processes, rather than just one.

In order to do estimation and verification of ocean-based carbon removal, both CDR suppliers and verifiers will need infrastructure to run earth systems models. The models themselves come from academia. ⟦C⟧worthy is a non-profit working on improving these models for CDR purposes.

Submarine is adapting these models to support the MRV process and building the data infrastructure needed to run them at scale for use by both CDR suppliers and verification services.

Let’s recap

Okay! That was a lot of information to swallow! Let’s do a quick recap.

To stabilize our climate, we need to dramatically reduce our emissions. But even with optimistic estimates for our rate of reduction, we’ll still need technology to remove CO2 directly from the atmosphere. To reach the scale of removal needed, we need to double the rate of global CO2 removal every 21 months for the next 27 years.

To make this happen, we need to accelerate a variety of market players and processes.

• Suppliers: Companies like Charm Industrial, Running Tide, Heirloom, Eion, Lithos, and Ebb Carbon have two main jobs: capturing carbon dioxide directly from the atmosphere, and storing it away forever in a way that doesn’t leak back into the atmosphere.
  • Permitting: To build industrial infrastructure, suppliers need permits. This process is long and painful right now. Companies like Blumen and Paces want to make this less painful.
  • Industry Partners: To gain access to land and permits, CDR suppliers partner with massive industries like agriculture, mining, steel manufacturing, and desalination plants. In exchange, the partners get decreased cost of waste removal, or co-benefits like improved soil quality.
• Capital: To fund operations, suppliers have a few key sources of capital:
  • Sale of carbon removal credits: Corporations and individuals voluntarily buy carbon removal credits. Companies like Watershed and Patch provide curated market access to enterprise customers. Companies like Wren and Commons provide access for individuals. Later, there may be regulatory requirements to buy these, or possibly direct procurement by governments.
  • Venture capital: Investment funds like Lowercarbon Capital, and MCJ Collective buy partial ownership of suppliers’ companies.
  • Grants: Organizations (mostly government agencies like the DOE or ARPA-E) provide money for research, which suppliers can apply for. To make this process easier, companies like Streamline and Pioneer are helping write and manage grant applications.
  • Loans: If suppliers can demonstrate that they will have money available to them down the road either from revenue or grant approvals, financial institutions like Enduring Planet will give them an advance on that capital.
• Trust: For the market to work, everyone involved must have consensus on how money is turning into permanently removed carbon dioxide.
  • Research: Academic, governmental, and industry research organizations publish papers on the fundamentals of a given CO2 removal method. These papers are peer reviewed and distributed by journals.
  • Life cycle analysis (LCA): Suppliers, usually in conjunction with a third party consultancy, publish a methodology showing the CO2 emissions and removal at every step in the process, demonstrating that the method creates net negative emissions.
  • Measurement, reporting, and verification (MRV): As suppliers run their process for removal in the real world, they carefully record measurements along the way. This can include things like the weight of bio-oil being injected, or the ocean salinity reading from a buoy where a fluid was released into the ocean. The LCA is then applied to determine the net removal from the actions of the company, which is entered into a registry so carbon credit buyers can see what their money was used for. Companies like Isometric and Submarine are working on improving these processes, alongside non-profit organizations like ⟦C⟧worthy.

Now then, what should you do with all this information?

Where to go from here

When asked what advice she’d offer to software engineers interested in contributing to the climate fight, here’s what the founder of Spark Climate Solutions said:

There are two overarching archetypes of paths I've seen here:

(a) find a software engineering job in a climate-related org or company, or

(b) go figure out where there are holes in the ecosystem, and use problem-solving, technical, and leadership skills to make a difference. https://t.co/uo8zqWnHn5

— Erika Reinhardt (@embrein) July 31, 2022

If you’re interested in path (a) and find yourself intrigued by carbon dioxide removal, many of the companies referenced in this article are hiring! Here are some job openings:

• Full Stack Developer @ Running Tide (Remote USA): Running Tide does CDR via marine biomass cultivation.
• Full Stack Software Engineer @ Ebb Carbon (San Carlos, California): Ebb Carbon does CDR via electrochemical ocean alkalinity enhancement.
• Lead Modeling and Simulation Engineer @ Heirloom (Brisbane, California): Heirloom does CDR via direct air capture.
• Technical Lead @ Undo (London, UK): Undo does CDR via enhanced rock weathering.
• Platform & SRE Engineering Roles @ Isometric (London, UK): Isometric is building a verification and multi-pathway registration service for carbon removal.
• Founding Engineer @ Streamline (San Francisco, California): Streamline helps climate companies apply for grants much faster.
• Founding Engineer @ Blumen (San Francisco, California): Blumen provides siting intelligence tools for geothermal, hydrogen, carbon sequestration, and minerals project developers.
• Multiple Software Engineering Positions @ Watershed (San Francisco, NYC, London, Remote): Watershed, among many other things, provides a curated marketplace for carbon removal purchases.
• Software Engineer @ Patch (San Francisco): Patch also does many things, but one is to provide a curated marketplace for carbon removal purchases.
• Full Stack Engineer, Climate @ Stripe. (Remote): Stripe helped catalyze the start of the permanent CDR market with early purchases.

If you’re interested in path (b), here’s the consistent advice I’ve gotten from people that have walked this path ⁹:

1. Talk to a bunch of people in a specific industry. If the problems are too diverse, narrow down the company category or job role until they become consistent.
2. Choose a specific problem. Ask people about the problem to understand if it’s a mild annoyance, or the source of a constant living hell. If you want to build a big venture capital-backed business, make sure solving this problem has a path to building a $1 billion company.
3. Propose a specific solution.
4. Try to rapidly figure out why your solution is, in fact, a bad idea. If you fail to find evidence it’s a bad idea, try to implement it. If you find strong evidence, return to step 2.

In the course of researching and writing this article, I’ve mostly been doing steps 1 and 2. Here are the problems I’m considering exploring for steps 3 and 4:

• MRV for ocean based CDR
• MRV for soil based CDR
• Removing duplication of effort for sensor data aggregation for new CDR companies

If you’re also interested in path (b), the good news is finding the frontier problems in the space doesn’t take that long, and people will be happy to see you when you get there:

A weird property of the frontier: finding the edge forces a realization that there’s very few people there, the others who’ve found it are tightly clustered and therefore quite happy to see you, and you all can’t help but ask “where is everybody?” in escalating confusion

— Ryan Orbuch (@orbuch) December 12, 2021

Further reading & acknowledgements

This article comes in a short lineage of folks building their understanding of CDR from the ground up, and sharing what they learned:

• “We Need To Take CO2 Out Of The Sky” by Ryan Orbuch. Ryan is now a Partner at Lowercarbon Capital leading their $350 million carbon removal fund. If you want to stay on top of the space, following him on Twitter is a good bet.
• “Scaling CDR” by Neil Hacker. Neil is now a Researcher at Isometric, helping to build their verification and registry service for carbon removal.

If you enjoyed reading this, here are some other resources you might like:

• For an alternate take on the current state of CDR, see The State of Carbon Dioxide Removal.
• For a crash course in CDR, see AirMiners’ free 6 week cohort-based course called Boot Up.
• For a video series on the ongoing evolution of industry, see OpenAir’s THIS IS CDR.
• For writing on how to shape the growing CDR market, see The Great Unwind.

Thanks to the folks at Lowercarbon Capital, MCJ Collective, Heirloom, Charm Industrial, Eion, Ebb Carbon, Isometric, Submarine, Vesta, Running Tide, Airminers, Blumen, Streamline, Avnos, Lillianah, and Lithos for sharing their knowledge and connections in the ecosystem.

Thanks to Ashley Zhang, Dhen Padilla, JN Fang, Jason Benn, Lyn Stoler, Max Krieger, Mishti Sharma, Neil Hacker, Temina Madon, and Ryan Gomba for reading drafts of this article and providing the invaluable feedback needed to shape it.

Special thanks to Owen Wang for being my exploration buddy throughout this endeavour, bouncing ideas around with me, and always reminding me to write down the questions I want to ask before each new meeting. I wouldn’t have had the energy to kickstart this without you.

If you’re interested in commissioning writing like this about a climate-related industry vertical, or about your own early-stage climate tech company, get in touch with me at jamie.lf.wong@gmail.com. If you’re building at the frontier of CDR (and especially if you’re working on MRV for mCDR), please reach out!

For the past few years for Christmas, I’ve asked people in my family for stories from their lives rather than a physical gift. Last year, my sister Emma, who works as a dentist, asked me for the same. This is the story I sent my sister, along with some illustrations by my very multitalented coworker Jessica Liu.

A lot of my day-to-day work as a software engineer involves debugging. We expect the system to do X but it does Y instead. Why?

One of my coworkers, Karl, really enjoys this process. He imagines the bug as a sneaky adversary, running around a maze, trying to evade him as he methodically searches through each turn. In the end, he always finds his hidden foe, wins a round of this game, and rests awaiting his next tour of the maze.

All of my most bizarre debugging adventures come from my time working at Figma. I’ve visited other offices, found bugs in Chrome, been bedevilled by a “Back to the Future” easter egg left in third party software, but the most outrageous lengths brought me into an hour long Lyft ride from San Francisco to Mountain View to meet a stranger who’d never heard of Figma on a Wednesday afternoon.

Debugging is often a process of making and testing hypotheses. Since you think the system should be doing one thing but it’s doing another, it must mean that one of your assumptions about how the system works is wrong.

For example, let’s say that your friend is trying to get into your apartment using your electronic keypad but it isn’t working. They might call you and tell you “Hey, I can’t get in.” You told them the code before, so they should be able to get in. Something you’re assuming is true clearly isn’t, so you ask them some questions to test your hypotheses.

“The code is 1499, is that what you’re entering?”
“Yeah, that’s what I entered. I’ll try it again right now… Nope. Still nothing.”

The first assumption you were testing was that they had the correct code. It turns out that assumption was correct, so that can’t be the problem. Time to test a different hypothesis.

“Do the buttons light up when you press them?”
“Yep. They light up. When I hit the second 9, the keypad turns red, blinks, and then goes dark.”

Okay, so the batteries aren’t dead. Your assumption that the keypad has power was right too. What could possibly be wrong?

“Uhh, are you sure you’re entering it in right?”
“Yes. Duh. I’m not an idiot. I press 1, then I press 4, then 9, then 9 again. It goes red. I turn the handle. Door doesn’t budge.”

“Hmm…”

“Are you sure you’re telling me the code right? My hands are freezing out here trying to enter this code, and your fireplace inside is just taunting me.”
“What fireplace? Wait. What address are you at?”
“907 Elk St, like you told me.”

You check the last text message you sent your friend. You did indeed send “907 Elk St”. You live at “907 Elm”.

Different bugs can require wildly different amounts of work to figure out.

One of the key things affecting how long it takes to diagnose and fix a bug is the duration of the feedback cycle. Figuring out the cause of a bug can take anywhere from one to hundreds of hypotheses tests to figure out the root cause of the issue. When it takes 100 cycles of “test hypothesis → form new hypothesis”, a feedback cycle of 10 seconds versus 10 minutes makes the difference between finding the cause in about 15 minutes and finding the cause in two days of full-time work.

All sorts of things can disrupt this feedback cycle.

Some bugs are tricky because they’re non-deterministic, which is a technobabble way of saying “even when you do the exact same things in the exact same order, you get different results.” Every once in a while, the laundry machine seems to swallow a sock, but you don’t know why. You think you might be able to stop it by always tying your socks together before putting them in the laundry machine. On the next load of laundry, you tie all your socks together in pairs and put them in. When you pull them out all your socks are still there! But did you fix the problem, or did you just get lucky this time? You can neither accept nor reject your hypothesis yet.

Another especially frustrating category of bug is nicknamed the “Heisenbug”, after the physicist Werner Heisenberg. Heisenberg was one of the first to assert that the mere act of observing a phenomenon may cause it to change its behaviour. In terms of software, what this means is that as soon as you start trying to change the code to have it tell you more about what it’s doing and what might be going wrong, the bug stops appearing! Imagine if your dentistry patient’s toothache went away when they opened their mouth even slightly.

But most of the gnarliest bugs I’ve had the imperative to work on at Figma suffer from a third kind of feedback loop problem: we can’t reproduce the problem, but some of our users can, and they hit it all the time.

This can happen for all sorts of reasons. They might have some Chrome extension installed that none of us do, and they didn’t realize that was an important part of the puzzle. Their network administrator might’ve blocked some communication channel that’s essential for Figma’s functioning over the internet. After diagnosing this kind of problem, we can usually either find a way of reproducing it by installing their software, or telling them to talk to their own coworkers to resolve the problem.

But the most frustrating version of the “no reproduction” issue is when something about their hardware is interacting poorly with our software. And the most fickle variety of hardware Figma deals with more directly than most other companies is graphics cards.

Graphics cards are the component in your computer that’s specially designed to efficiently perform the kinds of computations needed to update the millions of pixels on your screen at sixty frames per second.

I’m fond of computer graphics in part because the bugs can be so entertaining. Sometimes instead of text being drawn legibly on the screen, its upside-down and flickering madly. It’s got a special hall-of-mirrors kind of beauty to it.

These bugs are much less entertaining when they’re causing customers distress and we can’t reproduce the problems on our own hardware.

The first time I was involved in debugging a graphics card problem was in October, 2017. A designer at one of our early customers, a recruiting software company called Lever, reported that for some files, their entire canvas went black in Figma. We couldn’t reproduce, but the designer at Lever consistently could. Thankfully, Lever and Figma’s offices were a short bus ride away, so I went with our CTO Evan Wallace and one other coworker to investigate. They lent us their laptop (a MacBook Pro with an unusual graphics card inside), and I watched Evan test hypothesis after hypothesis about what could be going wrong on the machine using some special diagnostic software he’d written. After about 30 minutes, he was able to find a fix. We thanked our gracious hosts at Lever, and returned to our office, everyone quite happy.

The next time was in June, 2018. This time, a number of our users had written into us reporting that with a change in Figma’s graphics code, they started seeing patterns of visual noise appear in their design files. After chatting with many of the users through our support system, we were able to deduce that they all had the same series of graphics card. The GeForce GTX 10 series. This time, none of the affected users were anywhere in San Francisco, so we couldn’t just drop by their office. We could ask them to hop on a phone call, but it would be a lot to ask for them to sit with us for hours, typing what we asked them to into their computer word by word or giving us total remote access. Thankfully, I looked online and saw that Best Buy had laptops with this kind of graphics card in stock. So I took a car to the San Francisco Best Buy location, talked to one of the sales associates to confirm that the laptop I was looking had the card I needed, charged a few thousand dollars to the company card our CEO Dylan Field had lent me, and was on my merry way. A few hours later we’d figured out the cause and had a workaround to dodge the problem.

But by far the most outrageous instance of this was in May, 2019. We had just started rolling out changes I was spearheading to a core graphics component of Figma when a small portion of our user base wrote in to tell us something was wrong. Instead of the image taking up the full screen as they expected, it was about an eighth of the expected size, stuck up in the top left corner. After some careful analysis, we realized that they were all on Windows 7, and all using a specific version of a graphics card manufactured by Intel. Now knowing how this dance works, we looked to see if any of them were in San Francisco. Nope. One in Mumbai, one in Indonesia. We actually didn’t find any Figma users reporting this problem anywhere in the United States. Digging further, we discovered that this was a pretty old laptop, now out of vogue in America, but relatively common in Russia and South Asia. It was so old, in fact, that nobody sells them any more either.

Well, shit.

We knew that to diagnose this sanely, we needed to have the device in hand. My coworker Ryan mentioned a friend of his works in the device lab at Dropbox and they might be able to lend us a device. Large companies like Dropbox and Google keep large repositories of hardware configurations in their device lab to deal with exactly this contingency. But Figma was still under 100 people at the time, so we were far from having a fully stocked warehouse of obscure graphics cards.

While we waited to hear back from Ryan’s contact, the third teammate working with me on this graphics component, Lauren, turned to me and said “Hey Jamie, how many followers do you have on Twitter again?” “Just over 3000”, I responded. “Why don’t you see if any of them can help you?”

I thought there was no way that would work, but I had nothing to lose, so I sent out a tweet saying “Okay, I need to debug on a really specific old laptop configuration. I need a Window 7 or 8 laptop with an Intel HD Graphics card. If you have one and are in the SF bay area, I’ll buy you a meal or a drink if you can let me use it!”

Literally one minute later, a man named Scott replied “How long do you need it for?”

After confirming with him that he can reproduce the problem we’ve been seeing on this laptop, and that he’d be willing to lend me the laptop for a few weeks, I ask if he’d be willing to meet me in a coffeeshop somewhere convenient for him the next day. He agrees.

The next day I hopped in a Lyft from Figma’s office in San Francisco destined for a Starbucks in Mountain View. Once arrived at the Starbucks, I did the awkward “waiting for someone you’ve never met before” dance. Maybe you’ve done this before for a date where the profile picture was unclear, or maybe for a networking connection set up through a friend. After a few of the Starbucks patrons assured me that they were, indeed, not Scott, I sat down and waited.

Scott walked in a few minutes later. I stood up to greet him, and we made some idle chatter while I got us both drinks from the cashier. We settled into some tables nearby. As it turns out, the laptop Scott was going to lend me miraculously already had none of his own personal files on it. It was a spare laptop he kept around to drive custom karaoke events at different anime conventions. As long as I could return or replace the laptop before the event, he was happy to lend it to me. Scott works as a Technical Support Engineer at a big company, and just happened to follow me on Twitter after seeing my blog post about color. He’d never heard of Figma before, though was interested in it when I explained a bit more.

After talking for about half an hour, I thanked Scott profusely and returned to San Francisco in another hour long car ride. The day after, I was able to isolate the problem, create a workaround, and unblock our project. Reluctant to give up the laptop since it was the only way we had of debugging similar issues, I asked my manager for permission to expense sending a replacement laptop to Scott. The request was approved, and now Scott is the happy owner of a nicer, newer laptop, paid for by a company he had never heard of for deciding to help a person he had never met.

Scott is a great guy.

Thanks to Nikhil Thota, Spencer Chang, Andee Liao, and Lauren Budorick for providing feedback on drafts of this post, to Jessica Liu for the illustrations, and to Karl, Evan, Dylan, Ryan, and Lauren for being a part of this story!

You’re done your interviews, and you have a few Silicon Valley startup offers in-hand. The offers have a bunch of numbers and you’re not sure how to put them all together. You’re bombarded by terms like ISOs and NSOs and vesting period and cliff and strike price and fair market value that you don’t really understand, but owning 0.1% of the company sounds great in theory, so you go to sign on the dotted line.

Pause. Subtle differences here can have massive downstream consequences. If you’re unaware of some of these policies, you might end up paying $100,000+ more in taxes than you anticipated, or be forced to choose between staying at a job you don’t love and leaving potential millions on the table.

Back in 2012, I was at a conference called CUSEC. At the post-conference happy hour, an Engineering Manager at Twitter named Ian Chan took a bunch of us aside. He told us he was going to tell us something that could potentially save us a ton of money if we ever joined a startup. He proceeded to explain what an 83(b) election is, and he was right — it did and will save me a ton of money. I’d like to pass that knowledge on to you.

This post is going to focus on a narrow but common set of circumstances for employees at early startups (usually <$1B valuation). We’ll assume that the equity part of the offer is for incentive stock options (ISOs), rather than restricted stock units (RSUs), or actual stock if it’s a public company. Companies usually move from ISOs → RSUs → publicly tradable stock as the company grows. This guide is exclusively for ISOs. Since tax implications are important here, and taxes are specific to your locale, this post is also specifically references American taxes.

DISCLAIMER: I have zero professional experience in tax or law. I’m an early software engineer at Figma (company was ~20 people when I joined). Do not use this post as a tax calculator! Use this post a tool to ask better questions to well informed people about potential outcomes.

The standard deal

A standard deal for an early startup in the valley looks something like this:

• X shares
• Strike price of $Y
• 4 year vest, 1 year cliff

To make calculations easier for the rest of the post, I’m going to use the following semi-realistic numbers:

• 20,000 shares
• Strike price of $0.30
• 4 year vest, 1 year cliff

Let’s take a look at what each of these terms means.

20,000 shares: Your offer letter might say something like “Subject to the approval of the Board, you will be granted an option to purchase 20,000 shares”. This 20,000 is the number of shares from the company you’ll eventually be allowed to buy. This is similar to the # of shares you might buy of a public stock, except that you won’t be able to freely buy and sell these shares whenever you want to.

Strike price of $0.30: The strike price is the price you’ll pay per share to purchase. So if you want to buy 7,000 of the shares and your strike price is $0.30, then you’ll need to give the company 7,000⨉$0.30=$2,100. Crucially, as the value of the company increases (hopefully), this price does not increase. If the company increases in valuation from $1M to $50M between when you join and when you purchase, you still only pay $2,100, not $105,000. This increase in the value of the shares without an increase in the price to purchase is the whole point of stock options.

4 year vest, 1 year cliff: This is a standard “vesting schedule”. In simplest terms, this means if you stay for 4 years, you have the right to purchase the full number of shares, and if you stay for less than 1, you leave with nothing.

In more subtlety, it means that at the 1 year mark, you gain the right to purchase ¹⁄₄ of the total agreed upon share count (20,000/4=5,000). Every month after that, you’ll be able to purchase ¹⁄₄₈ of the remaining shares (20,000/48=416).

Once you hit each vesting date, you aren’t required to purchase. As long as you’re employed at the company, you have right to purchase any shares you’ve already vested.

It’s worth mentioning that many companies do an equity refresher program, where employees are granted additional equity after they’ve been at the company for a while. This is important to mitigate their total compensation decreasing dramatically once they’ve completed vesting their initial grant.

Time to purchase the stock

You’ve made it to your 1 year cliff! Congratulations! You’ve earned the right to purchase some of the company stock. Purchasing the stock is called called “exercising the option”. Since you have a four year vesting schedule, you now have the right to exercise ¹⁄₄ of your options. In the above offer, this allows you to purchase 20,000/4=5,000 shares. The purchase cost here is easy to calculate. It’s just the # of shares times the strike price. 5,000⨉$0.30 = $1,500.

But what about the taxes? Taxes can create an extremely employee-hostile situation here.

Let’s say that your company has been absolutely crushing it, and since you’ve joined, the company has increased its valuation 10x from $1M to $10M. The valuation that’s relevant here is called the “fair market value” or the “409A valuation”, which is sometimes different than the valuation in press releases about venture capital fund raises. Companies are required to get a new 409A valuation at least once every 12 months.

While you paid $0.30 per share, since the company has increased its valuation 10x, as far as the IRS is concerned, those exercised shares are worth 10x as well. So they’re worth $3.00 per share now. This sounds great until you realize that the IRS wants to tax you on the difference between those two amounts.

The exact tax details here are complicated, but the crucial thing is that you can get taxed by this thing called AMT (Alternative Minimum Tax) on the difference¹. Since we’re mostly interested in ballpark figures here, we’ll assume an AMT rate of 35% (28% Federal + 7% California State).

This means you’ll pay around (35%)⨉(5,000)⨉($3.00-$0.30)=$4,725 in taxes on top of the purchase price you pay to the company.

This brings your new total cash needed to $1,500+$4,725=$6,225. Notice how the taxes here are way more than the purchase price itself.

Now you might think “okay, well, that sucks, but I guess I can sell off some of my shares to cover the taxes on it”. Except you can’t. Even though you’ve earned the right to purchase the shares, there’s no legal market to sell the shares yet. So you can’t sell them. The shares are “illiquid”. So you’ve been heavily taxed on an asset that you can’t sell, and that might be worth $0 if the company goes bankrupt before you have any chance to sell.

Given that, a reasonable strategy might be to just hold onto the options until you can convert the shares into cash after you buy them. At that point, you will be able to sell shares to cover the taxes.

Time to quit

Given the scary tax treatment, you decide to wait to exercise. You just hold onto the options instead. You’ve now been at the company for 3 years, and you just got a really tempting offer from another company, so you’re considering quitting.

You might hope that you can just continue holding on to your stock options after you leave until the stock becomes liquid, but in many cases, you’ll be sorely disappointed.

Every stock option agreement will have what’s called a “post-termination exercise window”. This is the amount of time after you leave in which you’ll still be allowed to exercise your stock options. After this window closes, you forfeit all un-exercised stock options. The standard post-termination exercise window is only 90 days.

This leaves you with a choice between three crummy strategies:

1. Forfeit your equity
2. Pay heavy taxes on an illiquid asset that might ultimately be worth $0
3. Stay at the company, and turn down the other job offer

The same calculation as before applies to the exercise strategy, but just to lay it out, let’s assume the company has continued on a rocket trajectory and now has a fair market valuation of $100M, yielding a new fair market share value of $30.

You’ve now vested ³⁄₄ of your stock option agreement, allowing you to purchase (³⁄₄)⨉20,000=15,000 shares for 15,000⨉$0.30=$4,500. But the IRS now thinks these shares are worth $30.00 each, for a total value of 15,000⨉$30.00=$450,000.

Your AMT tax liability now becomes ($450,000-$4,500)⨉35% = $155,925. So the total amount of cash you need to put in to claim these shares becomes $4,500+$155,925=$160,420. This might be way more cash than you can afford to put in.

This situation is what people sometimes refer to as “golden handcuffs”. Given that their choices if they quit are “spend a ton of money on something very risky” or “walk away with nothing”, many people will simply stay until the company IPOs, even if they’d rather quit.

Avoiding the golden handcuffs

There are two main ways I know of which allows employees to avoid this situation: long post-termination exercise windows, and early exercise. These are both policies that companies may or may not offer to employees, so you should ask about them when considering your offer!

Long post-termination exercise window

Rather than 90 days, some companies offer much longer exercise windows. Figma has a 5 year window for all employees with a tenure of over 2.5 years. Coinbase has a 7-year post-termination exercise window for all employees with with a tenure of over 2 years. Pinterest had 7 before they IPO’d.

A 5 year post-termination exercise window means that you can leave the company and hold onto your options for up to 5 years after you leave. If any time during that 5 years the company’s shares become liquid, you’ll be able to exercise your options, buy the shares, and sell off some of them to cover your taxes. If, any time during that 5 years, the company goes bankrupt, you avoided spending any of your own money purchasing & paying taxes on an asset now worth $0.

There’s some subtlety here where the options change from ISOs to NSOs if you do this, which I won’t go into detail for. See the “ISO and NSO tax treatment” section of this blog post by Carta to learn more. At a high level it sounds like gains from ISOs are taxed as capital gains, and NSOs are taxed as income.

Early exercise

Another way to avoid the golden handcuffs problem is allow employees to purchase shares before the options are vested. Under an early exercise program, employees can purchase up to their entire grant immediately. If you leave the company before your vesting is complete, the company has the right to buy back the unvested shares at either the strike price or the current market value, whichever is lower.

So let’s say that you’ve just joined a company with the same 20,000 stock option, $0.30 strike price, 4 year vest, 1 year cliff stock option agreement, but that also has early exercise available.

You opt to exercise the entire agreement immediately. You pay 20,000⨉$0.30=$6,000 to the company immediately. To make this beneficial, you also need to file a form with the IRS called an “83(b) election”. This roughly tells the IRS that you want to be taxed on the equity now rather than when the stock options vest. You have to file your 83(b) within 30 days of exercising. Don’t mess this up.

What about the AMT? Well, you pay the difference between the current market value and the strike price. Assuming you exercise as soon as the company board sets your strike price, the difference should be zero, because the strike price is set at the fair market value! So if you get the timing right for this, you pay $0 in taxes at the point of exercise.

You still will get taxed on capital gains when you eventually sell the shares, but at least you can avoid getting taxed on the asset when it’s illiquid.

The downside of this is that you’re giving the company money to purchase an asset that might go down to zero at the point in time where you have the least information (the beginning). The upside is that you have no golden handcuffs! Whenever you choose to leave, you walk away with all of the equity you’ve vested to date without paying a dollar more.

Even if your company doesn’t have early exercise available when you join, it’s still worth discussing. Figma did not offer early exercise when I joined, but adopted the policy before the next 409A valuation, so I was still able to take advantage of it.

Even if you’ve had a new 409A valuation since your strike price was set, it’s possible that the difference between the fair market value and your strike price are still small enough that the tax liability is small, making early exercise still worthwhile.

Qualified Small Business Stock

There’s a more subtle benefit of exercising early if the company has gross assets under $50M: the stock may be considered Qualified Small Business Stock (QSBS). If it is considered QSBS, and if you’ve held the stock for 5 years or longer at the point of sale, then you’ll pay 0% capital gains to the IRS during that sale.

To see how much this matters, let’s say that the company IPOs 4 years into your tenure. It does well in the market, and 5 years after your start date, you decide you’d like to sell everything at the current price of $40. You have your entire stock option grant of 20,000 shares vested and exercised now, so the gross sales here will be 20,000⨉$40=$800,000, and the price paid was your strike price times the number of shares, so 20,000⨉$0.30=$6,000. For non-QSBS shares, you’ll be taxed federal capital gains, which is going to be at least 15% and possibly more. ($800,000-$6,000)⨉15%=$119,100. For QSBS shares, you pay no federal capital gains², so you save that entire $119,100! You’ll likely still owe state capital gains taxes, which can still be hefty, but you have to pay that regardless.

Questions you should be asking

If you’re weighing your different offers, here are some questions you should be asking about the equity to make sense of this. Each of these things is unlikely to be in the offer letter.

1. How long is the post-termination exercise window?
2. Is early exercise available? If not, why not? It’s possible the founders are unaware of this kind of policy existing, and would be happy to offer it. You can pitch it as a recruiting benefit to you and future potential employees. You can also push for this policy change even if you’ve already joined.
3. What is the strike price?
4. What do you think the total growth potential of the company is from current value? 10x? 100x?
5. What is the total number of shares outstanding? This is the “denominator” to think about when trying to understand what your number of shares means. This will let you calculate your % ownership of the company (ignoring dilution), which will help you gut check the maximum value of your shares assuming the company does really well.

Before I signed my Figma offer in 2016, Figma’s CEO Dylan Field said something akin to “do your financial planning assuming this equity ends up being worthless”. All of the examples in this post are incredibly optimistic, so keep in mind that everything could still hit zero.

Plan for the worst, hope for the best, and do the math.

Further reading

Carta, a company which manages startup equity, has similar guides with more details but less emphasis on how you can get screwed. They include descriptions of ISOs & NSOs and dilution which I skip in this post.

• Equity 101 part 1: Startup employee stock options
• Equity 101 Part 2: Stock option strike prices
• Equity 101 Part 3: How stock options are taxed
• What is a 409A valuation?

Thanks to Rudi Chen and Kevin Lau for feedback on drafts of this post, to Ian Chan for making me think about this in the first place, and to Julia Evans for writing a similar post about stock options which I’ve sent to a ton of people to in the past. Handcuffs Vector by Vecteezy.

I live alone at the moment, and the start of quarantine, was… rocky.

Here’s a jewel of an entry from my journal, from March 13:

Fuck, this week has been bad.

But, I think I’m kind of getting the hang of this. Here’s an entry from March 29:

I’m in a surprisingly good mood right now. I kind of have a headache, but aside from that I feel pretty good. Mood probably like an 8 / 10?

So I’d like to share a bit about what I think’s been helping bridge the gap between “Fuck, this week has been bad” to “8/10”. This isn’t a guide on how to be Your Best Self™️, or how to productivity hack your way into emerging from quarantine with a PhD and a six pack. This is a collection of ideas that I think have helped me reset myself to general okay-ness. First and foremost, this a post I’m writing as guidance for future-me who has temporarily forgotten all of this and is in need of it again.

Getting out of the pit

Earlier this year, I wrote about pulling myself out of the transient depths of despair aka the “Everything is Terrible” pit. I find myself looking for my list of tools for exiting the pit frequently. Early on in self-quarantine, I found myself in the same trap of feeling guilty about not being able to do more, but recognizing that I was in the pit helped me.

Here’s a continuation of my journal from the “Fuck, this week has been bad” entry:

I really want to be a force of positivity. I want to find ways to help my friends connect, and in a way that feels good to me too. To help pull us all out of this anxiety. But joining in a big group call just isn’t the thing that works for. I think that’s what being introverted really means to me. I need to connect with people 1:1, and I just need to do something with them that isn’t talk about COVID. But I don’t know how to get off the topic. I need to invest my time into something, anything that isn’t that. I definitely did notice myself slowly emerging from the pit after watching some Fullmetal Alchemist. I was able to clean my desk a little and put away dishes. And I guess I was able to pull myself out enough to be able to write this.

Okay, so I got out of the pit through recognizing I was there and taking the baby steps I needed to climb back into the cozy cabin. Next, let’s consider how to extend our stay in the cabin instead of tripping back into the pit every few steps.

Rebuilding my daily habits

Physical isolation from other people coinciding with the entire worldwide media cycle focusing (correctly) on a pandemic was a glorious kind of multi-targeted assassination of everyone’s emotional regulation toolkit. The exercise I get from going to the gym or dance classes? POOF. The ambient social interaction I get from having lunch with your coworkers? GONE. The warm glow of looking forward to parties or group trips? NOPE.

The removal of all of those things is out of my control. So, once I was back on steady ground momentarily, I started re-examining what was within my control to emotionally regulate. Starting with habits.

There’s a certain class of advice that’s:

1. Insulting cliche
2. Really well supported by science
3. Really easy to ignore, even after you’ve followed it in the past

For me, four of those standard nuggets of wisdom for emotional regulation are:

• Exercise regularly
• Meditate
• Sleep at the same time every day
• Talk to friends frequently

From the picture of my goal tracker below, you might be able to see some vague sort of pattern for how my mood correlated with my ability to complete my habits.

I track my habits on a piece of paper that I have clipped to my bedroom wall. Before you bombard me with your favourite apps for habit tracking, I like having this as a physical piece of paper. First, because I find scribbling in my habit completion way more satisfying that tapping a piece of glass, and second, because it’s an awful lot easier to avoid confronting an empty digital calendar once I start slipping than it is to avoid looking at my own bedroom wall forever.

In case you need to hear it, you can do this. You can build habits. You don’t have to choose any of the same habits for me, but hopefully you have some sense of what things are important for you to be happy, and hopefully some of them are still doable within the confines of your home while by yourself. List them out. Make your own grid on a piece of a paper. My pre-printed one filled up, and I don’t own a printer, so I made my own using a ruler and a pen:

Making a schedule

Having aspirational habits is all fine and good, but until you have a schedule that’s realistic about how much time each habit takes, it’ll be easy for them remain lofty aspirations.

My routine every day now looks like this:

7:00-7:30 wakeup, breakfast, watch anime
7:30-8:30 exercise while listening to podcasts
8:30-9:30 shower & meditate (I call a friend to meditate with)
10:00-6:00 work (except weekends)
6:00-7:00 dinner (I usually call a friend to talk while I cook)
7:00-10:00 different stuff every night
10:00-10:30 journal
10:30-11:30 get ready for bed & read
11:30-7:00 sleep

As it turns out, I almost always wake up 30-60 minutes later than I’m planning on and exercise a little less than I planned on and go to bed a little later than I was planning on, but I still get a fair chunk of this done.

I want to emphasize that these are the things I chose not because this is how I’m going to “win” at life, but because they seem to actually help me feel like me. If your routine looks more like this:

11:00-11:30 wakeup, breakfast, people watch through your window
11:30-12:00 doodle
12:00-1:00 order delivery from local restaurant, eat while watching 90 day fiancé**
1:00-6:00 free time! watch movies, game, read, doodle, TikTok
6:00-7:00 group call with friends eating leftovers from lunch
7:00-1:00 play Stardew Valley with friends
1:00-11:00 sleep

…and that makes you feel great, that’s great. Do that.

If you can understand what it is you need to be happy, and can schedule that so you don’t have to spend the energy every day to plan how to do them, then you’ll hopefully have more energy to actually do them.

If you have no idea what you need, experiment! Plan a schedule, do your best to stick to it for a week, and then re-evaluate at the end of the week.

Planning for the week

I noticed this really dumb pattern in my pre-shelter-in-place life. It would start when I wake up in the morning and say “Gee, I’d really like to have dinner with a friend tonight”. Then I look at the clock, discover I’m late to work, and bolt to catch the next subway to work. On the way, I open my phone to message friends, get distracted, and end up engaging in a work conversation on Slack. Then I’m at work, and every few hours I have the thought that I should really message a friend to grab dinner, but get pulled into something else then forget. Now it’s 6:00pm. I message a few friends, but they all unsurprisingly have plans for the meal that is now 30 minutes away. I grab takeout nearby, go home, and am a little bit sad.

Having dinner by myself when I wasn’t planning on it normally is a little demoralizing, but I get my social time in at work, so it’s no big deal. As far as I can tell, whether I had a good conversation with a friend is the single strongest indicator of whether I feel good at the end of the day during shelter in place, so it’s really important that I make this part of my day.

This is where batch planning really helps me out.

Every Sunday, I have a bunch of time blocked off to plan for the coming week.

The social planning part is how I avoid the last-minute-dinner-demoralization. I schedule out the coming week, and think about friends I haven’t caught up with in a while and set up calls over dinner. I’ll talk to them while I’m cooking most of the time.

I think this kind of batch planning has a few benefits over trying to schedule daily.

First, it gives me something to look forward to in the week! One of the most crushing things about Coronavirus for a lot of people has been the horrible lack of things to look forward to. So, make your own things to look forward to!

Second, it protects me against my own anti-social tendencies when I’m in the pit. When I feel like trash, I find it really hard to reach out. But if I batch plan this kind of stuff on Sunday and feel like a dumpster fire on a Tuesday, then Sunday-me’s got my back: friendship call at 6:00pm.

Lastly, I find that I just come up with more interesting ideas when I have explicit time set aside to plan stuff for the week. For example, I had this idea for a silly experiment where I asked a friend to send me a grocery list for a recipe without telling me what the recipe is. Tomorrow I’m going to call her and try to cook the recipe asking her only yes/no and numeric-answer questions and see how hilariously I screw it up. If I was trying to set up calls every day ad-hoc, I don’t think I would’ve stumbled on this kind of idea.

Old, nourishing media

Beyond the nothing-to-look-forward to conundrum, the early days of shelter-in-place were in part difficult for me because of how overwhelmingly difficult it was to think or talk about anything else. This is true not only for major news outlets, but also for nearly every single content producer I follow. Mark Rober, vlogbrothers, 3blue1brown, Smarter Every Day, Planet Money, etc. etc. It makes sense — once something is on everyone’s minds, it’s going to be difficult to produce any content that’s not that. The irony about feeling overwhelmed about this while writing content that is coronavirus topical isn’t lost on me.

Following the news to make well-informed decisions is a noble goal. But consuming a ton of media about crises tends to provide me a load of anxiety without actually compelling me into positive action. In the first few days, I was doing what Hank Green beautifully dubbed “The Anxious Scroll”:

When I’m doing the anxious scroll I feel as if I’m doing something useful, and I’m seeing the same 3-5 stories over and over again, so that they seem like 300 - 500 stories, and I want to know more and want to know what it’s gonna be like tomorrow, and I want to know what it’s gonna be like in 3 weeks, and I feel like I’m doing something that’s going to uncover that reality and uncover that truth and I’m just not.

So lately when I do consume media, I’ve been intentionally listening to old media. I finished rewatching Fullmetal Alchemist: Brotherhood. I’ve been listening to the old Radiolab mini-series G about intelligence and the theft of Albert Einstein’s brain. I’ve been rewatching Avatar: The Last Airbender with some friends.

I’m not advocating for sticking your head in the sand and pretending that the world is business as usual. Stay informed enough to make sensible decisions both for yourself and your community, but stay away from the anxious scroll.

Laughing at myself

I think there’ve been spans of my life where I take myself and my problems altogether too seriously. While I’m a strong proponent of confronting your problems rather than minimizing them by trying to brush them off as inconsequential compared to the problems of others (this helps nobody), I also think there’s a beautiful light in being able to find humour in dark places.

To that effect, I’ll leave you with a little story.

In a group chat I’m in, a few of my friends started exchanging videos of them with their partners or quarantine buddies attempting to do some TikTok challenges, like this:

@blake.wood

Seesaw challenge completed ✅ ##fyp ##xyzcba ##seesaw ##MagicBoots @lorna.bensted

♬ Lalala - Y2K & bbno$

In watching my friends excitedly exchange their attempt videos, I felt a little sad.

I don’t have a romantic partner to do this with at the moment, and I don’t have a quarantine buddy to make an attempt.

What I do have, however, is Photoshop.

Hang in there, even if the only physical support you have around is yourself.

Thanks to Andee Liao, Spencer Chang, and Lilly Shi for reading drafts and encouraging me to post this.

Okay, Jamie, you’re in the pit again. It’s far from your first stay here, but let’s not get too comfortable. There’s a distinct lack of roses here to stop and smell. Mostly just thorns.

I know it’s not easy. Your body is lead. Your mind is in a slow, murky swirl. Your emotional energy level is at about a 3%.

Let’s go through your move list and take stock.

Move #277: Reach out to a close friend to talk? Sounds promising in net energy, but looks like the activation energy is more like a 10% if you have to be the one to reach out. Rough. Looks like we’ve got to go way down the activation energy list here to gain some ground.

Move #86: Clean your room maybe? Ah, closer, but that sounds more like a 5% kind of task. There’s just so much stuff to clean up. Seems like a lot of work.

Move #435: Read social media? Only requires 1% activation energy! Great! Start reading, I’ll check back in 30.

Hey, Jami — Oh… I see. You’re still in bed, just… just staring at the ceiling. Cooooool.

Let’s get your energy reading. 2%?! You went down a percent? Huh. Guess I didn’t look at the net energy prediction for that social media binge too carefully. Not great.

I’m not even sure what moves exist down here. Obviously the social media move is still there, but that didn’t go so well, so let’s explore some other options.

Let’s see, let’s see. Damn, most of the standard options are just totally out of range right now. Move #451: Read a book, needs 20% activation. Move #555: Go to the gym, needs 15%.

How ab — oh good you’re listening to music now. Wait. Wait, what music is this. Oh.

Move #393: Listen to music that reminds you of your exes. Requirement: 1%. Well, I guess that’s a marginal improvement from staring at the ceiling.

It’s… it’s working? I wouldn’t say you look happy exactly, but the meter says 5%!

Holy shit folks, he’s left the bed! 👏

Dang you’ve got your combo rolling now. #86: Clean your room. Done. Yield: +2%. #222: Shower. Done. +1%. #16: Water your plant. +1%. Boom. You’ve got some momentum. Meter reads 9%.

Your general vibe reads melancholy but I can see a dim life-is-worth-living glint in your eyes if I look past the dark circles.

After you make your morning smoothie and rewatch an episode of Fullmetal Alchemist (Brotherhood, of course), you’re at a cool 11% energy. I think I even saw you smile briefly. Time for you to plan ahead a little. Let’s revisit #277: Reach out to a friend to talk.

You can do it Jamie, just pick some people and fire off a couple of messages. “You free for dinner tonight?” That’s it. These are people that’ve known you for years. Nothing to fear here.

Ahhh you’re back in bed, huh? Glued to Reddit. 7% and… falling. Friends all already have dinner plans and you kicked yourself for not making plans earlier? Yeah, we’ve seen this play-by-play before.

Come on dude, let’s go get some sunlight and some groceries. Maybe try #44: Listen to some comedy. I think Ali Wong’s Baby Cobra is on Spotify now.

Okay, okay, now we’re talking! A little bit of laughter, a little bit of Vitamin D glow, and ingredients for a healthy lunch. Let’s get it! This 15% energy reading seems like a great foundation for things to come. You’re a little fragile still, but your options are opening up.

Alright, I’m going to drop that narrative device before it gets too old. The above, second-person narrative fell out of a framework that congealed in my mind today.

In almost all situations, I want to gain emotional energy and spread that to everyone around me. Every action I take has some effect on my energy level. But the effect is also going to depend on what energy level I’m at right now. And some actions are just not reasonably available to me at certain energy levels.

While thinking about specific percentages has a certain charm, spending a lot of time trying to dial in whether something has an expected value of +4% or +6% isn’t really the point of the framework to me. So let’s talk about three ranges of energy I find myself in: the pit, the cozy cabin, and the dance floor.

The entire story opening this post takes place in the “Everything is Terrible” pit. When I’m in this state, it’s a struggle to do anything. When I’m this low energy, there is zero creative juice flowing. Focus is stolen by mind fog and rumination. Staying in bed feels both inescapable and shameful.

So when I’m in this state, what I need is not something challenging. I need easy wins. I need reminders that I am, in fact, a reasonably competent human capable of enjoyment and basic task completion. For me, this means consuming familiar media like The Anthropocene Reviewed, Key & Peele, and Ali Wong, or just doing really basic maintenance tasks around my apartment like doing laundry.

Once I pull myself out of the pit, I’ve hopefully entered the “Life is Pretty Good” cozy cabin. All things considered, this place is nice. I’ve got some hot chocolate spiked with some Bailey’s in my mug, the fire is crackling, some jazzy Christmas music is playing. This chair is SO. SOFT. I could spend a good long time here without getting sick of it.

From the cozy cabin, whole categories of actions that felt insurmountably challenging start to seem pretty straightforward. I might schedule dinners with friends for the rest of the week, or head to the gym, or take a dance class. And those actions will typically prolong my cabin stay.

But hopefully those actions will also start to unblock a further set of actions: the risky ones that pave the path to the dance floor. If the cabin is the place where I feel at peace, the dance floor is the place where I feel alive.

Almost by definition, the actions leading to the “I AM THE SHIT” dance floor have to be risky. When I’m on this metaphorical dance floor, I’m in a state of supreme confidence that only comes from overcoming some sort of fear.

And by definition of “risk”, sometimes the gamble won’t pay off. I ask someone out and get ghosted. I go to a house party and don’t connect. I take the challenging dance class and totally blank when it’s my turn to perform. When this happens, if I’m lucky, I just get bounced to the cabin. If it goes particularly poorly, I might find myself back in the pit.

I like this model because it acknowledges that the things I need when I’m low energy are very different from what I need when I’m high energy. From the pit, I frequently feel frustrated that the higher energy activities feel out of reach, and I’m hoping that reminding myself of this model will be calming.

By thinking about this model explicitly, I was forced me to write down what I think I need in each of these places. I’m hoping that next time I’m in the pit, I can look at this and use it to help me crawl out a little bit faster.

This is where I’d normally try to put some pithy quotation or closing thought tying the whole thing together, but I’ve got three other posts that are in 90%-completion hell, and I kind of need an easy win today. So I’ll just leave you with a question: what does your pit → cabin → dance floor transition plan look like? Send me a picture, I wanna see!

Thanks to Owen Wang for helping pull me out of the pit tonight and for coaching me how to better season my rice & beans over the phone. I wouldn’t have had the emotional or caloric energy to write this without you.

When we’re little kids, we naturally channel our emotional state as a stream of consciousness. We cry and scream at a pin drop: when we’re hungry, when we’ve been away from a parent’s embrace for too long, when we’re sleepy, when there’s poop in our diapers, when we have to share our favourite toy but we just don’t wanna. We’re skilled in allowing our emotions to flow, but totally hopeless at guiding the path of the emotional expression. We’re untethered emotional firehoses. This is an ineffective route of expression.

As we mature, we learn that unfettered, uncontrolled expression of emotion is ineffective. We see that people pull away from us when we yell at them. We see that throwing tantrums in public embarrasses the people we’re with. So for many people, like me, we learn to suppress. We turn the emotions inwards, convince ourselves that they’re useless, and believe we can restrain them with the right kind of patience and resolve. In doing so, we avoid the kind of external damage caused by screaming at people, or flying into a violent rage, or breaking our hands punching lockers. But we also compress what should be fleeting moods into something deeper and darker. It might not explode, but that suppression and deeply held hurt pervades our daily experience, taking the form of shame, disconnection, and fear. This is the route of suppression.

So as we mature still, we must learn to find channels to allow that emotion to run its course, but do so constructively. I still struggle a lot with this.

We release through appreciating art which reflects our pain. We release through creating art which reflects our pain. We release through resolving conflict by returning to it with honesty about how we feel and what triggered us to feel that way. This is an effective route of expression.

When I’m upset with my girlfriend, I could take the ineffective route of expression: “You’re 15 minutes late! You always do this. You don’t care about me”.

Or, I can take the route of suppression: proceed with dinner, trying to pretend I was never angry. Telling myself that I had no right to be angry.

Or, I can take a more effective route of expression: “I felt angry when I saw that I had been waiting 15 minutes past the time we agreed upon for you to come to dinner. Next time, would you please either try harder to be on time, or try to let me know how late you’re going to be?” Beginning the dialog this way invites empathy, rather than defensiveness. Once this happens, I might discover the anger relents to reveal the more vulnerable fear that underlies it. And with the empathetic support of a partner, that fear too might be allowed to run its course. And from there, security and playfulness arise.

Going to a therapist by myself, I learned a bit about how to release and integrate past pain by allowing long held, deeply held emotion to flow.

Going to see a couples therapist with my girlfriend, we’ve learned a bit about how to release and integrate past pain we’ve inflicted upon ourselves and each other by allowing each other to see our emotions flow. By providing dedicated space for emotionally charged discussions about our relationship, we see the pain felt by one another. We learn to understand it and empathize with it. And we learn to have emotionally honest conversations outside of the therapy sessions. I think that since we’ve started going to therapy, we don’t have conflicts less frequently, but we reach resolution to conflict much faster when it arises. And each time after we reach resolution, we’re able to return to warmth and play, uninhibited by distance from unspoken conflict.

So if you’re curious, and if the resources are available to you as an individual or as a pair, go to a therapist and feel those feels.

Hello little evil
In the depths of my chest
You’ve taken shelter too long

How easy for you to hide in the shadows
Of an unexamined space
Living in your home
And giving nothing back

I know you mean no harm
You are not malicious
You simply are

How unfortunate for you
When a stray rock shatters a tinted window
Of your shrouded home
And you are brought to light

Blinded
You thrash, You sob, You scream
As you are evicted
And it echos in your host

I wonder when the next rock
Will happen upon your now departed home
And reveal your brethren

April 25, 2014

Hey you,

It’s me: you. I made it. I have my iron ring¹. I guess you aren’t worried about failing to get here, but here looks a little different than you probably hoped it would. You’re not going to work at Google, and even more surprisingly, you won’t want to. You’ll discover that in your own time, so let me calm some of your fears and instill in you a few more.

In high school, you were the techie among your friends, and that was cool. It made you different and interesting and helped you create cool things to show your friends. In university, you’ll be surrounded by people with similar skill sets and interests. Your classmates will catch up to you. You will panic a little when this happens. Calm down. You’re not competing with them. Respect them, and learn from them, and never compare your accomplishments with theirs.

Make friends within your program. They’ll help you get through material that doesn’t click in your head as quickly as you’re used to. You won’t get integration, electromagnetics, or feedback control the way you got most things in high school. They’ll also end up being your teachers and your roommates. They’ll keep you motivated to do all your work, even though they won’t motivate you to go to class.

Oh yeah – you’re not going to go to all your classes. Thinking less of the people who skip in first year is going to make you feel pretty stupid when you almost completely stop going to class in upper years. You’ll stop going initially because you learn the material faster by reading textbooks in the library, but eventually you’ll just focus on learning through assignments and cram when you need to. This is maybe not ideal, but it’ll free you to get more exercise.

Make friends outside your program. As odd as this might sound to you now, you’re going to get sick of being surrounded by people in tech 24 / 7. Talk to a diverse group of people. Socializing in university will become a conscious thing. You won’t make friends during lectures like you did in class during high school.

Never stop playing badminton. Join the executive team as early as possible. This will be one of the most transformative things you do in university. Going out for dinner with them will make you feel more socially confident. The friends you make through badminton club will make you feel like part of a community in a very different way than your class does, and a few of them will be critical to meeting the incredible girl I’m with today.

Get more comfortable talking to strangers. Whenever you go meet a group of people, you’re not going to get along with all of them. That’s okay. Calling for pizza delivery over the phone should not be a source of social anxiety. Traveling will help with this immensely.

Open up to people. To get people to open to you, sometimes you have to go first. Find people where this comes naturally for you. When you do find them, don’t hold back, and make sure to let them know you care. Talking to them in groups of people will be very different than talking to them 1-on-1. This is true of everyone - especially your family. Do what you can to stay in touch with your family and get to know them. There’s a lot more going on in their lives than you think. Schedule Skype dates and dinners with family and friends.

Be humble. People will learn about your achievements in their own way, and if they don’t that’s okay. Being approachable and fun to talk to is far more important than deference. Stop wearing tech shirts and hoodies all the time. It makes you look simultaneously aesthetically oblivious and less approachable.

Don’t stop reading. Read non-fiction, science fiction, and some of the classics. Books will alter your perception of the world, lead to some incredibly engaging conversations, and strengthen future friendships. They’ll also give you something interesting to think about during rough periods in your life. You’ll extract more meaning reading a book for 20 minutes than 2 hours of skimming Hacker News. You should still look up this Hacker News thing though, it’s pretty handy.

Experiment with what makes you productive. You’ll think at first that to-do lists work for you, but they really don’t. You work better on pretty much everything in the morning, before you check your email, go on Facebook, or read Hacker News. This is especially true when you have to read a lot of notes, or when you want to write anything at all. Let yourself take naps while studying.

Keep a journal. It will let you reminisce, maintain your confidence, and organize your thoughts in difficult times. No matter when you started, you’ll wish you started earlier. You’ll find the act of writing itself cathartic, even when you do have people to tell your worries to.

Face your problems head on, but use the support you have. Even if your problems seem insignificant compared to others’, that doesn’t give you an excuse to ignore them. If you’ve been thinking the same thought on and off for years, it’s not going to go away by itself.

Oh, and when you go ATV riding, if you feel like you’re shifting around a lot, you’re doing it wrong. Hug the chassis with your thighs. Trust me on this one, you’ll save your ankle a world of hurt.

— Jamie Wong, University of Waterloo Software 2014 Alumnus

An edited version of this post was written for Mozilla Hacks: Cross-language Performance Profile Exploration with speedscope.

For the past 9 months, I’ve been working on speedscope: a fast, interactive, web-based viewer for large performance profiles. You can use it live at www.speedscope.app, and read the code on GitHub at jlfwong/speedscope.

It’s inspired by the performance panel of Chrome developer tools and by Brendan Gregg’s FlameGraphs. If you’ve never heard of flamegraphs before or have heard of them but never understood how to read them, the guide “Using flamegraphs” from rbspy’s documentation is wonderful.

The goal of speedscope is to provide a 60fps way of interactively exploring large performance profiles from a variety of sources. It runs totally in-browser, and does not send any profiling data to any servers. Because it runs in-browser, it should work in Chrome and Firefox on Mac, Windows, and Linux (though I’ve only actually tested on Mac 😬).

You can use it offline via the command-line if you install it via npm.

npm install -g speedscope

Or you can just download a zipfile from the GitHub releases page and open the contained index.html directly if you don’t want to deal with node or npm.

Currently supported import sources include:

• JavaScript
  • Importing from Chrome
  • Importing from Firefox
  • Importing from Node.js
• Ruby
  • Importing from stackprof
  • Importing from rbspy
• Native code
  • Importing from Instruments.app (macOS)
  • Importing from perf (linux)

speedscope also provides two ways of importing from custom sources. Check out Importing from custom sources for the details. If you’d like to see support for an additional format, please submit an issue, and ideally contribute that support!

speedscope is MIT licensed.

What can it do?

speedscope is broken down into three primary views: Time Order, Left Heavy, and Sandwich.

🕰Time Order

In the “Time Order” view (the default), call stacks are ordered left-to-right in the same order as they occurred in the input file, which is usually going to be the chronological order they were recorded in. This view is most helpful for understanding the behavior of an application over time, e.g. “first the data is fetched from the database, then the data is prepared for serialization, then the data is serialized to JSON”.

The horizontal axis represents the “weight” of each stack (most commonly CPU time), and the vertical axis shows you the stack active at the time of the sample. If you click on one of the frames, you’ll be able to see summary statistics about it.

⬅️Left Heavy

In the “Left Heavy” view, identical stacks are grouped together, regardless of whether they were recorded sequentially. Then, the stacks are sorted so that the heaviest stack for each parent is on the left – hence “left heavy”. This view is useful for understanding where all the time is going in situations where there are hundreds or thousands of function calls interleaved between other call stacks.

🥪 Sandwich

The “Sandwich” view is a table view in which you can find a list of all functions and their associated times. You can sort by self time or total time.

It’s called the “Sandwich” view because if you select one of the rows in the table, you can see flamegraphs for all the callers and callees of the selected row.

Who’s using it?

I haven’t done too much to advertise speedscope, but it’s already in use in at least three companies!

Figma

We’ve been using speedscope at Figma for all kinds of performance analysis for our user interface design tool. This is a bit of a cheat because I work at Figma, but it spread from initially a tool that only I used to a tool that’s in use by most people on the teams that work in C++ or TypeScript for performance analysis. One particular superpower that speedscope gives us is the ability to get users to record performance profiles in Chrome containing minified names, then remap the stack traces back to development symbols using the .symbols file generated by emscripten. We’ve been using it for months now to help us optimize loading time in particular, and regularly use it to load 100MB+ profiles exported from Chrome.

Benchling

Engineers at Benchling have been using it thanks to Alan Pierce using it to implement a really easy-to-use profiler for their backend Python code.

We recently worked on a project to make backend profiling easier and more accessible, and speedscope has been an important part of that. Now anyone in the company (not just engineers) can take a profile anytime they see something slow in any production instance, post a speedscope link to Slack, and anyone on engineering can help out in the analysis and share insights or ideas for improvements. As one recent example, speedscope recently made it easy to see that a slow operation spent almost all of its time building and running regexes, and switching to Python’s string.find when possible sped the operation up by about 75x.

Uber

Rudro Samanta and teammates have been using speedscope at Uber to visualize profiles of iOS applications captured via Xcode & Instruments.

Instruments is great at finding a specific hotspot function across a long running trace of our app (e.g. Font loading – UIFont fontWithName is slow). With speedscope we’re able to find that during this specific transition from one screen to the other, this specific icon load appears to take up more time than we expected and maybe we should pre-load that one large icon or something.

Architecture

While working on speedscope, I tried to stay generally light on dependencies, opting to implement a lot of stuff from scratch rather than trying to glue a bunch of third party libraries together. This serves two major purposes:

1. Small download size. At time of writing, the total compressed download size of https://www.speedscope.app/ is 131KB, and only the first 63KB of that is needed before the page becomes usefully interactive. This fully loads in just under 2.5s over a Fast 3G connection.
2. Re-inventing the wheel teaches me stuff. This might not be the best strategy for getting product in front of users as fast as possible, but I certainly learn a lot more doing this. As part of speedscope, I implemented a basic linear algebra math library, pan/zoom mechanics, a least-recently-used cache, and an efficient list view, among other things. I feel much better about how I spent the time designing & implementing those than I would have if I spent that time learning the specific APIs of other people’s implementations.

That said, there were a few dependencies that were gnarly enough that I didn’t want to roll my own.

• Preact for declaratively defining UI views. I generally really enjoy working with React, but wanted a smaller library in terms of byte size since one of my goals was to make speedscope load really fast even from an empty cache on a slow connection.
• Aphrodite for authoring CSS. This deals with all the auto-prefixing for me, avoids classes of bugs around CSS prioritization conflicts, and works nicely with TypeScript, which I’ve really thoroughly enjoyed using for this project. I have a bias here vs. other CSS solutions since I co-authored the first version of Aphrodite while at Khan Academy, so I understand how it works ☺️.
• Redux and preact-redux for global in-memory state management. I honestly started trying to implement my own global state management solution for speedscope since redux felt like a lot of unnecessarily complexity for what I was trying to do, but I repeatedly ran into bugs or complications that led me back to redux. I ultimately left with a greater appreciation for redux.
• Pako for extracting zlib encoded data as part of importing trace files from Instruments. I probably would’ve learned a lot implementing my own decompressor, but it seems like something pretty easy to get wrong, and it’s a pretty small dependency.
• graphics.ts is a WebGL abstraction layer that I ported from https://github.com/evanw/sky. speedscope uses WebGL for flamegraph rendering, which is how it ends up being so fast even for massive profiles. Before that I was using regl, which was delightful for the most part. I ultimately needed to switch away from it because it uses eval, which would preclude speedscope from being usable in environments with a strict Content-Security-Policy header. See https://github.com/regl-project/regl/issues/491#issuecomment-414146457.

As a quick note about build and test infrastructure, speedscope uses Jest for writing tests, Travis CI for continuous integration, Coveralls for test code coverage reporting, Parcel for source code → build artifact transformation, and prettier for automatic code formatting.

What’s next?

Working on speedscope has provided me a bounty of blogpost material. Working on it led me to writing both “Color: From Hexcodes to Eyeballs” and “Reverse Engineering Intruments’ File Format”, and there’s a lot more material I could extract from this project to write about.

In terms of spreading speedscope to wherever it can be useful, there are two major paths: increasing the number of formats that speedscope can import, and making tighter integrations into speedscope to make the import process seamless. For instance, I have an open PR on stackprof to make it the default visualization tool that ships with the profiler.

In terms of increasing the capabilities of speedscope itself, one path forward would be to support a wider array of visualizations to e.g. be able to visualize network requests & CPU time on the same time axis, just like Chrome developer tools does. In an extreme version of this, you could imagine importing complex traces from distributed systems and exploring them interactively.

Buuuuuut I think I’m probably going to take a long break before trying to do any of that (though I’m still happy to do code review if someone else wants to!) Let me explain why.

At the beginning of the year, I finished reading 7 Habits of Highly Effective People, and wrote down my personal mission statement broken down by the different roles I care about: friend, son, brother, student, teacher, etc. Part of what led me to focusing on speedscope for the last 9 months was this part of my mission statement:

As an engineer, I will work to build maintainable, understandable systems and tools which aid people in solving problems of essential complexity in the world.

I feel confident that speedscope does contribute to that goal in a small way. But I think that I’ve let too much of my time drift towards this. While I like working on speedscope, it’s a mostly solitary endeavor, and even if I was collaborating with other people, I already spend ~40 hours doing satisfying engineering work at Figma. I think I’d like to reclaim more of my time to pursuing creative, collaborative endeavors with other people where most of that collaboration doesn’t happen through a sheet of glass. So that’s likely where most of my time is going to go for the rest of this year. This stuff comes in waves, so perhaps I’ll return to sinking a ton more time into this some time in the future.

Part of what I like about side projects and having my own blog is not needing to keep working on stuff when I don’t really feel like it, and not needing to pretend that the work is my lifeblood. While I’m proud of a lot of the technical stuff I’ve written on my blog, I’m ultimately most proud of my post about my experience with depression, because it expresses a distinctly human perspective that I rarely take the time to express.

Welp, this “product launch” post took a weirdly personal turn near the end. If you got this far, here’s a music video from one of my favorite artists 😃.

Have you ever wondered how applications store their data? Plenty of file formats like MP3 and JPG are standardized and well documented, but what about custom, proprietary file formats? What do you do when you want to extract data that you know is in a file somewhere, and there are no APIs to extract it?

Over the last few months, I’ve been building a performance visualization tool called speedscope. It can import CPU profile formats from a variety of sources, like Chrome, Firefox, and Brendan Gregg’s stackcollapse format.

At Figma, I work in a C++ codebase that cross-compiles to asm.js and WebAssembly to run in the browser. Occasionally, however, it’s helpful to be able to profile the native build we use for development and debugging. The tool of choice to do that on OS X is Instruments. If we can extract the right information from the files Instruments outputs, then we can construct flamecharts to help us build intuition for what’s happening while our code is executing.

Up until this point, all of the formats I’ve been importing into speedscope have been either plaintext or JSON, which lends them to easier analysis. Instruments’ .trace file format, by contrast, is a complex, multi-encoding format which seems to use several hand-rolled binary formats.

This was my first foray into complex binary file reverse engineering, and I’d like to share my process for doing it, hopefully teaching you about some tools along the way.

Disclaimer: I got stuck many times trying to understand the file format. For the sake of brevity, what’s presented here is a much smoother process than I really went. If you get stuck trying to do something similar, don’t be discouraged!

A brief introduction to sampling profilers

Before we dig into the file format, it will be helpful to understand what kind of data we need to extract. We’re trying to import a CPU time profile, which helps us answer the question “where is all the time going in my program?” There are many different ways to analyze runtime performance of a program, but one of the most common is to use a sampling profiler.

While the program being analyzed is running, a sampling profiler will periodically ask the running program “Hey! What are you doing RIGHT NOW?”. The program will respond with its current call stack is (or call stacks, in the case of a multithreaded program), then the profiler will record that call stack along with the current timestamp. A manual way of doing this if you don’t have a profiler is to just repeatedly pause the program in a debugger and look at the call stack.

Instruments’ Time Profiler is a sampling profiler.

After you record a time profile in Instruments, you can see list of samples with their timestamps and associated call stacks.

This is exactly the information we want to extract: timestamps, and call stacks.

Exploring with file and tree

If you’d like to follow along with these steps, you can find my test file here: simple-time-profile.trace, which is a profile from Instruments 8.3.3. This is a time profile of a simple program I made specifically for analysis without any complex threading or multi-process behaviour: simple.cpp.

A good first step when trying to analyze any file is to use the unix file program.

file will try to guess the type of a file by looking at its bytes. Here are some examples:

$ file favicon-16x16.png
favicon-16x16.png: PNG image data, 16 x 16, 8-bit colormap, non-interlaced
$ file favicon.ico
favicon.ico: MS Windows icon resource - 3 icons, 48x48, 256-colors
$ file README.md
README.md: UTF-8 Unicode English text, with very long lines
$ file /Applications/Google\ Chrome.app/Contents/MacOS/Google\ Chrome
/Applications/Google Chrome.app/Contents/MacOS/Google Chrome: Mach-O 64-bit executable x86_64

So let’s see what file has to say about our .trace file.

$ file simple-time-profile.trace
simple-time-profile.trace: directory

Interesting! So the Instruments .trace file format isn’t a single file, but a directory.

macOS has a concept of a bundle, which is effectively a directory that can act like a file. This allows many different file formats to be packaged together into a single entity. Other file formats like Java’s .jar and Microsoft Office’s .docx. accomplish similar goals by grouping many different file formats together in a zip compressed archive (they’re literally just zip archives with different file extensions).

With that in mind, let’s take a look at the directory structure using the tree command, installed on my Mac via brew install tree.

$ tree -L 4 simple-time-profile.trace
simple-time-profile.trace
├── Trace1.run
│   ├── RunIssues.storedata
│   ├── RunIssues.storedata-shm
│   └── RunIssues.storedata-wal
├── corespace
│   ├── MANIFEST.plist
│   ├── currentRun
│   │   └── core
│   │       ├── extensions
│   │       ├── stores
│   │       └── uniquing
│   └── run1
│       └── core
│           ├── core-config
│           ├── extensions
│           ├── stores
│           ├── table-manager
│           └── uniquing
├── form.template
├── instrument_data
│   └── 20202640-0B46-4698-ADAD-DF54B3ABE816
│       └── run_data
│           └── 1.run.zip
├── open.creq
└── shared_data
    └── 1.run

…okay then! There’s a lot going on in here, and it’s not clear where we should be looking for the data we’re interested in.

Finding strings with grep

Strings tend to be the easiest kind of data to find. In this case, we expect to find the function names of the program somewhere in the profile. Here’s the main function of the program we profiled:

int main(int argc, char* argv[]) {
  while (true) {
    alpha();
    beta();
    gamma();
    delta();
  }
  return 0;
}

If we’re lucky, we’ll find the string gamma somewhere in plaintext in the .trace bundle. If the data were compressed, we might not be so lucky.

$ grep -r gamma simple-time-profile.trace
Binary file simple-time-profile.trace/form.template matches

Cool, so form.template contains the string gamma in it somewhere. Let’s see what kind of file this is.

$ file simple-time-profile.trace/form.template
simple-time-profile.trace/form.template: Apple binary property list

So what’s this Apple binary property list thing?

Interpreting the plist with plutil

From a Google search, I found an article about converting binary plists, which references a tool called plutil for analyzing and manipulating the contents of binary plists. plutil -p seems especially promising as a way of printing plists in a human readable format.

$ plutil -p simple-time-profile.trace/form.template
{
  "$version" => 100000
  "$objects" => [
    0 => "$null"
    1 => "rsrc://Template - samplertemplate"
    2 => {
      "NSString" => <CFKeyedArchiverUID ...>{value = 3}
      "NSDelegate" => <CFKeyedArchiverUID ...>{value = 0}
      "NSAttributes" => <CFKeyedArchiverUID ...>{value = 5}
      "$class" => <CFKeyedArchiverUID ...>{value = 11}
    }
    ...(many more entries here, excluded for brvity)
  ]
  "$archiver" => "NSKeyedArchiver"
  "$top" => {
    "com.apple.xray.owner.template" => <CFKeyedArchiverUID ...>{value = 12}
    "com.apple.xray.instrument.command" => <CFKeyedArchiverUID ...>{value = 234}
    "$1" => <CFKeyedArchiverUID ...>{value = 163}
    "cliTargetDevice" => <CFKeyedArchiverUID ...>{value = 0}
    "com.apple.xray.owner.template.description" => <CFKeyedArchiverUID ...>{value = 2}
    "$2" => <CFKeyedArchiverUID ...>{value = 164}
    "com.apple.xray.owner.template.version" => 2.1
    "com.apple.xray.owner.template.iconURL" => <CFKeyedArchiverUID ...>{value = 1}
    "$0" => <CFKeyedArchiverUID ...>{value = 141}
    "com.apple.xray.run.data" => <CFKeyedArchiverUID ...>{value = 247}
  }
}

I wasn’t familiar with many Mac APIs, so the best I could do was just Google search some of the terms in here. CFKeyedArchiverUID shows up a lot here, and that sounds related to NSKeyedArchiver.

A Google search tells me that NSKeyedArchiver is an Apple-provided API for serialization and deserialization of object graphs into files. If we can figure out how to reconstruct the object graph that was serialized into this, this might be instrumental in extracting the data we need!

A convenient way to explore Cocoa APIs is inside of an XCode Playground. Inside of a playground, I was able to construct an NSKeyedUnarchiver and start to pull data out of it, but I quickly ran into problems:

In particular, we have this error:

cannot decode object of class (XRInstrumentControlState) for key (NS.objects);
the class may be defined in source code or a library that is not linked

Unsurprisingly, in order to decode the objects stored within a keyed archive, you need to have access to the classes that were used to encode them. In this case, we don’t have the class XRInstrumentControlState, so the archiver has no idea how to decode it!

We could probably work around this limitation by subclassing NSKeyedUnarchiver and overriding the method which decides which class to decode to based on the class name, but I ultimately want to be able to read these files via JavaScript in the browser where I won’t have access to the Cocoa APIs. Given that, it would be helpful to understand how the serialization format works more directly.

To be able to extract data from this file, we’ll both need to be able to do the same thing as plutil -p is doing above, and also do the same thing as a NSKeyedUnarchiver would be doing in reconstructing the object graph from the plist file.

Making a binary plist parser

Thankfully, parsing binary plists is a problem that many others have encountered in the past. Here are some binary plists parsers in a variety of languages:

• JavaScript: node-bplist-parser
• Java: BinaryPropertyListParser.java
• Python: bplist.py
• Ruby: rbBinaryCFPropertyList.rb

Ultimately, I ended up making minor modifications to a binary plist parser that we use at Figma for Sketch import, which you can now find in the speedscope repository in instruments.ts.

Reconstructing the object graph

Also helpfully, other people have done analysis on how NSKeyedArchiver serializes its data into a property list. This blogpost by mac4n6, for example, explores an example of how an object graph can be reconstructed from the property list. It ends up being a relatively straightforward process of replacing numerical IDs with their corresponding entries in the $object lookup table.

After these replacements are completed, objects will have a property indicating what class they were serialized from. Many common datatypes have consistent serialization formats that we can use to construct a useful representation of the original object.

This too ends up being a task that surprisingly many people have been interested in solving, and have also kindly release source code to solve:

• Python: ccl_bplist.py and bpylist.py
• JavaScript: msUnarchiver.js

You can see examples of these in patternMatchObjectiveC in the speedscope source code.

Handling custom datatypes

There are datatypes in simple-time-profile.trace/form.template, however, that are specific to Instruments. When we’re trying to reconstruct an object from an NSKeyedArchive, we’re given a $classname variable. If we collect all the Instruments-specific classnames and print them out, we’re left with this:

[
  "XRRecordingOptions",
  "XRContext",
  "XRAnalysisCoreDetailNode",
  "XRAnalysisCoreTableQuery",
  "XRMainWindowUIState",
  "XRInstrumentControlState",
  "XRRunListData",
  "XRIntKeyedDictionary",
  "PFTPersistentSymbols",
  "XRArchitecture",
  "PFTSymbolData",
  "PFTOwnerData",
  "XRCore",
  "XRThread",
  "XRBacktraceTypeAdapter"
]

Stepping back, what we’re trying to figure out here is where the function names and file locations are stored within this file. From surveying the list of Instruments specific classes above, PFTSymbolData seems like a good candidate to contain this information.

A google search of PFTSymbolData yields this github page showing a reverse-engineered header file from XCode!

@interface PFTSymbolData : NSObject <NSCoding, CommonSymbol, XRUIStackFrame, NSCopying>
{
    NSString *sourcePath;
    struct XRLineNumData *addressData;
    int numAddresses;
    int addressesCapacity;
    BOOL _missingSymbolName;
    struct _CSRange symbolRange;
    unsigned int fTypeFlags;
    PFTOwnerData *ownerData;
    NSMutableArray *inlinedInstances;
    NSString *symbolName;
}

These headers were extracted using class-dump. This was pretty lucky — it just so happens that someone has dumped all of the headers in XCode and put it up in a GitHub repository.

Using the header as a reference and inspecting the data, I was able to reconstruct a semantically useful representation of PFTSymbolData

You can see the relevant code in readInstrumentsKeyedArchive.

Now we have the symbol table, but we still need the list of samples!

Finding the list of samples with find and du

I was hoping that all of the information I was interested in would be in a single file within the .trace bundle, but it turns out we aren’t so lucky.

The next thing I’m looking for is the list of samples collected during instrumentation. Each sample contains a timestamp, so I expect them to be stored as a table of numbers. But I wasn’t even sure what numbers I should be looking for, because I’m not sure how the timestamps are stored. The timestamps could be stored as absolute values since unix epoch, or could be stored as relative to the previous sample, or relative to the start of the profile, and could be stored as floating point values or integers, and those integers might be big endian or small endian.

Overall, I wasn’t really sure how to find data when I didn’t know any values that would definitely be in the data table, so I had a different idea for an approach. I recorded a longer profile, then went looking for big files! I figured that as profiles got longer, the data storing the list of samples should get bigger.

To find potential files of interest, I ran the following unix pipeline:

$ find . -type f | xargs du | sort -n | tail -n 10
152     ./corespace/run1/core/stores/indexed-store-9/spindex.0
168     ./corespace/run1/core/stores/indexed-store-12/spindex.0
296     ./corespace/run1/core/stores/indexed-store-3/bulkstore
600     ./form.template
808     ./corespace/run1/core/stores/indexed-store-9/bulkstore
1064    ./corespace/run1/core/stores/indexed-store-12/bulkstore
2048    ./corespace/run1/core/uniquing/arrayUniquer/integeruniquer.data
2048    ./corespace/run1/core/uniquing/typedArrayUniquer/integeruniquer.data
20480   ./corespace/currentRun/core/uniquing/arrayUniquer/integeruniquer.data
20480   ./corespace/currentRun/core/uniquing/typedArrayUniquer/integeruniquer.data

Let’s break down this pipeline.

• find . -type f finds all files in the current directory, printing them one per line (find man page)
• xargs du runs du to find file size using the list piped to it as arguments using xargs. We could alternatively do du $(find . -type f). (xargs man page, du man page)
• sort -n numerically sorts the results in ascending order (sort man page`)
• tail -n 10 takes the last 10 lines of output (tail man page)

We can extend this command to tell us the file type of each of these files:

$ find . -type f | xargs du | sort -n | tail -n 10 | cut -f2 | xargs file
./corespace/run1/core/stores/indexed-store-9/spindex.0:                     data
./corespace/run1/core/stores/indexed-store-12/spindex.0:                    data
./corespace/run1/core/stores/indexed-store-3/bulkstore:                     data
./form.template:                                                            Apple binary property list
./corespace/run1/core/stores/indexed-store-9/bulkstore:                     data
./corespace/run1/core/stores/indexed-store-12/bulkstore:                    data
./corespace/run1/core/uniquing/arrayUniquer/integeruniquer.data:            data
./corespace/run1/core/uniquing/typedArrayUniquer/integeruniquer.data:       data
./corespace/currentRun/core/uniquing/arrayUniquer/integeruniquer.data:      data
./corespace/currentRun/core/uniquing/typedArrayUniquer/integeruniquer.data: data

The cut command can be used to extract columns of data from a plaintext table. In this case cut -f2 selects only the second column of the data. We then run the file command on each resulting file. (cut man page)

A file type of data isn’t very informative, so we’ll have to start examining the binary contents to figure out the format ourselves.

Exploring binary file contents with xxd

xxd is a tool for taking a “hex dump” of a binary file (xxd man page). A hex dump of a binary file is a representation of a file displaying each byte of the file as a hexadecimal pair.

$ echo "hello" | xxd
00000000: 6865 6c6c 6f0a                           hello.

The output here shows the offset (00000000:), the hex represenation of the bytes in the file (6865 6c6c 6f0a) and the corresponding ASCII interpretation of those bytes (hello.), with . being used in place of unprintable characters. The . in this case corresponds to the byte 0a which in turn corresponds to the ASCII \n character emitted by echo.

Here’s another example using printf to emit 3 bytes with no printable representations.

$ printf "\1\2\3" | xxd
00000000: 0102 03                                  ...

Let’s use this to explore the biggest file we found.

$ xxd corespace/currentRun/core/uniquing/typedArrayUniquer/integeruniquer.data | head -n 10
00000000: 6745 2301 7e33 0a00 0100 0000 0000 0000  gE#.~3..........
00000010: 0100 0000 0000 0000 0000 0000 0000 0000  ................
00000020: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000030: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000040: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000050: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000060: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000070: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000080: 0000 0000 0000 0000 0000 0000 0000 0000  ................
00000090: 0000 0000 0000 0000 0000 0000 0000 0000  ................

Hmm. It seems like there’s a lot of data in this file that’s all zero’d out. The sample list can’t possibly be all zeros, so we’d rather just look for bits that aren’t zero’d out. The -a flag of xxd can be of help in this situation.

       -a | -autoskip
              toggle autoskip: A single '*' replaces nul-lines.  Default off.


$ xxd -a ./corespace/currentRun/core/uniquing/typedArrayUniquer/integeruniquer.data
00000000: 6745 2301 7e33 0a00 0100 0000 0000 0000  gE#.~3..........
00000010: 0100 0000 0000 0000 0000 0000 0000 0000  ................
00000020: 0000 0000 0000 0000 0000 0000 0000 0000  ................
*
009ffff0: 0000 0000 0000 0000 0000 0000 0000 0000  ................

Welp. It doesn’t seem like this file actually contains much useful data. Let’s re-sort our list of files, this time sorting by the number of non-null lines.

$ for f in $(find . -type f); do echo "$(xxd -a $f | wc -l) $f"; done | sort -n | tail -n 10
    1337 ./corespace/currentRun/core/extensions/com.apple.dt.instruments.ktrace.dtac/knowledge-rules-0.clp
    1337 ./corespace/run1/core/extensions/com.apple.dt.instruments.ktrace.dtac/knowledge-rules-0.clp
    2232 ./corespace/currentRun/core/extensions/com.apple.dt.instruments.poi.dtac/binding-rules.clp
    2232 ./corespace/run1/core/extensions/com.apple.dt.instruments.poi.dtac/binding-rules.clp
    2391 ./corespace/run1/core/uniquing/arrayUniquer/integeruniquer.data
    2524 ./corespace/run1/core/stores/indexed-store-12/spindex.0
    2736 ./corespace/run1/core/stores/indexed-store-9/spindex.0
    5148 ./corespace/run1/core/stores/indexed-store-9/bulkstore
    6793 ./corespace/run1/core/stores/indexed-store-12/bulkstore
   18091 ./form.template

We already know what form.template is, so we’ll start with the second largest. If we look at indexed-store-12/bulkestore, it looks like there might be some useful data in there, starting at offset 0x1000.

$ xxd -a ./corespace/run1/core/stores/indexed-store-12/bulkstore | head -n 20
00000000: 0a0a 3412 0300 0000 2800 0000 0010 0000  ..4.....(.......
00000010: 2100 0000 0040 0800 0040 0000 0000 0000  !....@...@......
00000020: 0000 0000 0000 0000 0000 0000 0000 0000  ................
*
00001000: 796c 8f2b 0000 0400 0000 0000 0006 0000  yl.+............
00001010: 0004 0000 0040 420f 0000 0000 00fe 0000  .....@B.........
00001020: 00a4 9ddc 2b00 0004 0000 0000 0000 0200  ....+...........
00001030: 0000 0400 0000 4042 0f00 0000 0000 0001  ......@B........
00001040: 0000 5d7f 0a2c 0000 0400 0000 0000 0002  ..]..,..........
00001050: 0000 0004 0000 0040 420f 0000 0000 0002  .......@B.......
00001060: 0100 0039 fb19 2c00 0004 0000 0000 0000  ...9..,.........
00001070: 0000 0000 0400 0000 4042 0f00 0000 0000  ........@B......
00001080: 0401 0000 336b 292c 0000 0400 0000 0000  ....3k),........
00001090: 0000 0000 0004 0000 0040 420f 0000 0000  .........@B.....
000010a0: 0005 0100 0026 e538 2c00 0004 0000 0000  .....&.8,.......
000010b0: 0000 0000 0000 0400 0000 4042 0f00 0000  ..........@B....
000010c0: 0000 0701 0000 4c5a 482c 0000 0400 0000  ......LZH,......
000010d0: 0000 0000 0000 0004 0000 0040 420f 0000  ...........@B...
000010e0: 0000 0009 0100 0040 ce57 2c00 0004 0000  .......@.W,.....
000010f0: 0000 0000 0000 0000 0400 0000 4042 0f00  ............@B..

The @B in the right column, while not obviously semantically meaningful, seems to repeat at a regular interval. Maybe if can figure out that regular interval, we’ll be able to guess what the structure of the data is. We can try guessing different intervals by using the -c argument of xxd, and try changing the byte grouping using the -g argument.

      -c cols | -cols cols
              format <cols> octets per line. Default 16 (-i: 12, -ps: 30, -b: 6).

     -g bytes | -groupsize bytes
              separate the output of every <bytes> bytes (two hex characters or eight bit-digits each) by a whitespace.  Specify -g 0 to suppress grouping.

We seem to get alignment between the repeated values when we group data into chunks of 33:

$ xxd -a -c 33 -g0 ./corespace/run1/core/stores/indexed-store-12/bulkstore | cut -d' ' -f2 | head -n 10
0a0a34120300000028000000001000002100000000400800004000000000000000
000000000000000000000000000000000000000000000000000000000000000000
*
00000000796c8f2b000004000000000000060000000400000040420f0000000000
fe000000a49ddc2b000004000000000000020000000400000040420f0000000000
000100005d7f0a2c000004000000000000020000000400000040420f0000000000
0201000039fb192c000004000000000000000000000400000040420f0000000000
04010000336b292c000004000000000000000000000400000040420f0000000000
0501000026e5382c000004000000000000000000000400000040420f0000000000
070100004c5a482c000004000000000000000000000400000040420f0000000000

Sweet! This suggests that this file format uses 33 bytes per entry, and hopefully each of those entries corresponds to one sample in the profile.

Looking around in the directory which contains this bulkstore, we find a helpful sounding file called schema.xml:

$ cat ./corespace/run1/core/stores/indexed-store-12/schema.xml
<schema name="time-profile" topology="XRT50_C22_TypeID">
    <column engineeringType="XRSampleTimestampTypeID" engineeringName="Sample Time" mnemonic="time" topologyField="XRTraceRelativeTimestampFieldID"></column>
    <column engineeringType="XRThreadTypeID" engineeringName="Thread" mnemonic="thread" topologyField="XRCategory1FieldID"></column>
    <column engineeringType="XRProcessTypeID" engineeringName="Process" mnemonic="process"></column>
    <column engineeringType="XRCPUCoreTypeID" engineeringName="Core" mnemonic="core"></column>
    <column engineeringType="XRThreadStateTypeID" engineeringName="State" mnemonic="thread-state"></column>
    <column engineeringType="XRTimeSampleWeightTypeID" engineeringName="Weight" mnemonic="weight"></column>
    <column engineeringType="XRBacktraceTypeID" engineeringName="Stack" mnemonic="stack"></column>
</schema>

Alright, this is looking pretty good! XRSampleTimestampTypeID and XRBacktraceTypeID seem particularly relevant.

The next step is to figure out how these 33 byte entries map onto the fields in schema.xml.

Guessing binary formats with Synalyze It!

So far in this exploration, all of the tools I’ve been using come standard on most unix installations, and all are free and open source. While I certainly could have figured this out end-to-end using only tools in that category, my friend Peter Sobot introduced me to a tool that made this process much easier.

Synalyze It! is a hex editor and binary analysis tool for OS X. There’s a Windows & Linux version called Hexinator. These tools let you make guesses about the structure of file formats (e.g. “I think this file is a list of structs, each of which is 20 bytes, where the first 4 bytes are an unsigned int, and the last 16 bytes are a fixed-length ascii string”), then parse the file based on that guess and display in both a colorized view of the hex dump and in an expandable tree view. This let me guess-and-check several hypotheses about what the structure of the file.

Eventually I was able to guess the length and offsets of the fields I was interested in. Synalyze It! helps you visually parse the information by setting colors for different fields. Here, I’ve set the sample timestamp to be green, and the backtrace ID to be red.

From looking at the values of the sample time and comparing them with what Instruments was displaying, I was able to infer that the values represented the number of nanoseconds since the profile started as a six byte unsigned integer. I was able to verify this by editing the binary file and then re-opening it in instruments.

Sweet! So that answers the question of where the sample information is stored, and we know how to interpret the timestamp data. But we still don’t quite know how to turn the backtrace ID into a stack trace.

To try to find the stacks, we can see if the memory addresses identified as part of the symbol table show up anywhere outside of the form.template binary plist.

Finding binary sequences using python

Here’s the same symbol data from earlier.

So let’s see if we can find one of these addresses referenced somewhere else in the .trace bundle. We’ll look for the third address in that list, 4536213276.

As a first attempt, it’s possible that the number is written as a string.

$ grep -a -R -l '4536213276' .

No results. Well, that was kind of a long shot. Let’s try the more plausible idea of searching for a binary encoding of this number.

There are two standard ways of encoding multi-byte integers into a byte stream. One is called “little endian” and the other is called “big endian”. In little endian, you place the least significant byte first. In big endian, you place the most significant byte first. Using python’s struct standard library, we can see what each of these representations look like.

The number is too big to fit in a 32 bit integer, so it’s probably a 64 bit integer, which would make sense if it’s a memory address that has to support 64 bit addresses.

$ python -c 'import struct, sys; sys.stdout.write(struct.pack(">Q", 4536213276));' | xxd
00000000: 0000 0001 0e61 1f1c                      .....a..
$ python -c 'import struct, sys; sys.stdout.write(struct.pack("<Q", 4536213276));' | xxd
00000000: 1c1f 610e 0100 0000                      ..a.....

>Q instructs struct.pack to encode the number as a big endian 64 bit unsigned integer, and <Q corresponds to a little endian 64 bit unsigned integer.

If you split up the bytes, you can see it’s the same bytes in both encodings, just in reverse order:

Now we can use a little python program to search for files with the value we care about.

$ cat search.py
import os, glob, struct

addr = 4536213276
little = struct.pack('<Q', addr)
big = struct.pack('>Q', addr)

for (dirpath, dirnames, filenames) in os.walk('.'):
  for f in filenames:
    path = os.path.join(dirpath, f)
    contents = open(path).read()
    if little in contents:
      print 'Found little in %s' % path
    elif big in contents:
      print 'Found big in %s' % path
$ python search.py
Found big in ./form.template
Found little in ./corespace/run1/core/uniquing/arrayUniquer/integeruniquer.data

Sweet! The value is in two places: one little endian, one big endian. The form.template one we already knew about — that’s where we found this address in the first place. The second location in integeruniquer.data is one we haven’t explored. It also was one of the files we found when searching for files with large amounts of non-zero data in them.

After fumbling around in this file with Synalyze It! for a while, I discovered that file is aptly named, and contains arrays of integers packed as a 32 bit length followed by a list of 64 bit integers.

So integeruniquer.data contains an array of arrays of 64 bit integers. Neat! It seems like each 64 bit int is either a memory address or an index into the array of arrays. This was the last piece of the puzzle we need to parse the profiles.

Putting it all together

So overall, the final process looks like this:

1. Find the list of samples by finding a bulkstore file adjacent to a schema.xml which contains the string <schema name="time-profile".
2. Extract a list of (timestamp, backtraceID) tuples from the bulkstore
3. Using the backtraceID as an index into the array represented by arrayUniquer/integeruniquer.data, convert the list of (timestamp, backtraceID) tuples into a list of (timestamp, address[]) tuples
4. Parse the form.template binary plist and extract the symbol data from PFTSymbolData from the resulting NSKeyedArchive. Convert this into a mapping from address to (function name, file path) pairs.
5. Using the address → (function name, file path) mapping in conjunction with the (timestamp address[]) tuple list, construct a list of (timestamp, (function name, file path)[]) tuples. This is the final information needed to construct a flamegraph!

Phew! That was a lot of digging for what ultimately ends up being a relatively straightforward data extraction. You can find the implementation in importFromInstrumentsTrace in the source for speedscope on GitHub.

If you do get the chance to give speedscope a try, please tweet @jlfwong and let me know what you think 🙂.

Thanks to Peter Sobot, Ryan Kaplan, and Rudi Chen for providing feedback on the draft of this post.

This post is also available in Russian: Цвет: от шестнадцатеричных кодов до глаза, and Japanese: 色：ヘキサコードから眼球まで.

Why do we perceive background-color: #9B51E0 as this particular purple?

This is one of those questions where I thought I’d known the answer for a long time, but as I inspected my understanding, I realized there were pretty significant gaps.

Through an exploration of electromagnetic radiation, optical biology, colorimetry, and display hardware, I hope to start filling in some of these gaps. If you want to skip ahead, here’s the lay of the land we’ll be covering:

Otherwise, let’s start with the physics.

Electromagnetic radiation

Radio waves, microwaves, infrared, visible light, ultraviolet, x-rays, and gamma rays are all forms of electromagnetic radiation. While these things all go by different names, these names really only label different ranges of wavelengths within the electromagnetic spectrum.

The smallest unit of electromagnetic radiation is a photon. The energy contained within a photon is proportional to the frequency of its corresponding wave, with high energy photons corresponding with high frequency waves.

To really understand color, we need to first understand radiation. Let’s take a closer look at the radiation of an incandescent light bulb.

We might want to know how much energy the bulb is radiating. The radiant flux ($$\Phi_e$$) of an object is the total energy emitted per second, and is measured in Watts. The radiant flux of a 100W incandescent lightbulb is about 80W, with the remaining 20W being converted directly to non-radiated heat.

If we want to know how much of that energy comes from each wavelength, we can look at the spectral flux. The spectral flux ($$\Phi_{e,\lambda}$$) of an object is radiant flux per unit wavelength, and is typically measured in Watts/nanometer.

If we were to graph the spectral flux of our incandescent lightbulb as a function of wavelength, it might look something like this:

The area under this curve will give the radiant flux. As an equation, $$\Phi_e = \int_0^\infty \Phi_{e,\lambda}(\lambda) d\lambda$$. In this case, the area under the curve will be about 80W.

Now you might’ve heard from eco-friendly campaigns that incandescent lightbulbs are brutally inefficient, and might be thinking to yourself, “well, 80% doesn’t seem so bad”.

And it’s true — an incandescent lightbulb is a pretty efficient way to convert electricity into radiation. Unfortunately, it’s a terrible way to convert electricity into human visible radiation.

Visible light

Visible light is the wavelength range of electromagnetic radiation from $$\lambda = 380\text{nm}$$ to $$\lambda = 750\text{nm}$$. On our graph of an incandescent bulb, that’s the shaded region below.

Okay, so the energy radiated within the visible spectrum is $$8.7W$$ for an efficiency of $$8.7\%$$. That seems pretty awful. But it gets worse.

To understand why, let’s consider why visible light is, well, visible.

Human perceived brightness

Just as we saw that an incandescent light bulb doesn’t radiate equally at all wavelengths, our eyes aren’t equally sensitive to radiation at all wavelengths. If we measure a human eye’s sensitivity to every wavelength, we get a luminosity function. The standard luminosity function, $$\bar y(\lambda)$$ looks like this:

The bounds of this luminosity function define the range of visible light. Anything outside this range isn’t visible because, well, our eyes aren’t sensitive to it!

This curve also shows that our eyes are much more sensitive to radiation at 550nm than they are to radiation at either 650nm or 450nm.

Other animals have eyes that are sensitive to a different range of wavelengths, and therefore different luminosity functions. Birds can see ultraviolet radiation in the range between $$\lambda=300\text{nm}$$ to $$\lambda=400\text{nm}$$, so if scholarly birds had defined the electromagnetic spectrum, that would’ve been part of the “visible light” range for them!

By multiplying the graph of spectral flux with the luminosity function $$\bar y(\lambda)$$, we get a function which describes the contributions to human perceived brightness for each wavelength emitted by a light source.

This is the spectral luminous flux ($$\Phi_{v,\lambda}$$). To acknowledge this is about human perception rather than objective power, the luminous flux is typically measured in lumens rather than Watts using a conversion ratio of 683.002 lm/W.

The luminous flux ($$\Phi_v$$) of a light source is the total human perceived power of the light.

Just as we calculated the radiant flux by taking the area under the spectral flux curve, we can find the luminous flux by taking the area under the spectral luminous flux curve, with a constant conversion from perceived watts to lumens:

So the luminous flux of our 100W incandescent lightbulb is a measly 2.4W or 1600lm! The bulb has a luminous efficiency of 2.4%, a far cry from the 80% efficiency converting electricity into radiation.

Perhaps if we had a light source that concentrated its emission into the visible range, we’d be able to get more efficient lighting. Let’s compare the spectra of incandescent, fluorescent, and LED bulbs:

And indeed, we can see that far less of the radiation in fluorescent or led bulbs is wasted on wavelengths that humans can’t see. Where incandescent bulbs might have an efficiency of 1-3%, fluorescent bulbs can be around 10% efficient, and LED bulbs can achieve up to 20% efficiency!

Enough about brightness, let’s return to the focus of this post: color!

Quantifying color

How might we identify a given color? If I have a lemon in front of me, how can I tell you over the phone what color it is? I might tell you “the lemon is yellow”, but which yellow? How would you precisely identify each of these yellows?

Armed with the knowledge that color is humans’ interpretation of electromagnetic radiation, we might be tempted to define color mathematically via spectral flux. Any human visible color will be some weighted combination of the monochromatic (single wavelength) colors. Monochromatic colors are also known as spectral colors.

For any given object, we can measure its emission (or reflectance) spectrum, and use that to precisely identify a color. If we can reproduce the spectrum, we can certainly reproduce the color!

The sunlight reflected from a point on a lemon might have a reflectance spectrum that looks like this:

Note: the power and spectral distribution of radiation that reaches your eye is going to depend upon the power & emission spectrum of the light source, the distance of the light source from the illuminated object, the size and shape of the object, the absorption spectrum of the object, and your distance from the object. That’s a lot to think about, so let’s focus just on what happens when that light hits your eye. Let’s also disregard units for now to focus on concepts.

When energy with this spectral distribution hits our eyes, we perceive it as “yellow”. Now let’s say I take a photo of the lemon and upload it to my computer. Next, I carefully adjust the colors on my screen until a particular point of the on-screen lemon is imperceptibly different from the color of the actual lemon in my actual hand.

If you were to measure the spectral power distribution coming off of the screen, what would you expect the distribution to look like? You might reasonably expect it to look similar to the reflectance spectrum of the lemon above. But it would actually look something like this:

Two different spectral power distributions that look the same to human observers are called metamers.

To understand how this is possible, let’s take a look at the biology of the eye.

Optical biology

Our perception of light is the responsibility of specialized cells in our eyes called “rods” and “cones”. Rods are predominately important in low-light settings, and don’t play much role in color vision, so we’ll focus on the cones.

Humans typically have 3 kinds of cones. Having three different kinds of cones makes humans “trichromats”. There is, however, at least one confirmed case of a tetrochromat human! Other animals have even more cone categories. Mantis shrimp have sixteen different kinds of cones.

Each kind of cone is labelled by the range of wavelengths of light they are excited by. The standard labelling is “S”, “M”, and “L” (short, medium, long).

These three curves indicate how sensitive the corresponding cone is to each wavelength. The highest point on each curve is called the “peak wavelength”, indicating the wavelength of radiation that the cone is most sensitive to.

Let’s see how our cones process the light bouncing off the lemon in my hand.

By looking at the normalized areas under the curves, we can see how much the radiation reflected from the real lemon excites each of cones. In this case, the normalized excitations of the S, M, and L cones are 0.02, 0.12, and 0.16 respectively. Now let’s repeat the process for the on-screen lemon.

Despite having totally different radiation spectra reaching the eye, the cone excitations are the same (S=0.02, M=0.12, L=0.16). That’s why the point on the real lemon and the point on the digital lemon look the same to us!

Our 3 sets of cones reduce any spectral flux curve $$\Phi_e(\lambda)$$ down to a triplet of three numbers $$(S, M, L)$$, and every distinct $$(S, M, L)$$ triplet will be a distinct color! This is pretty convenient, because (0.02, 0.12, 0.16) is much easier to communicate than a complicated continuous function. For the mathematically inclined, our eyes are doing a dimensional reduction from an infinite dimensional space into 3 dimensions, which is a pretty damn cool thing to be able to do subconsciously.

This $$(S, M, L)$$ triplet is, in fact, our first example of a color space.

Color spaces

Color spaces allow us to define with numeric precision what color we’re talking about. In the previous section, we saw that a specific yellow could be represented as (0.02, 0.12, 0.16) in the SML color space, which is more commonly known as the LMS color space.

Since this color space is describing the stimulation of cones, by definition any human visible color can be represented by positive LMS coordinates (excluding the extremely rare tetrachromat humans, who would need four coordinates instead of three).

But, alas, this color space has some unhelpful properties.

For one, not all triplet values (also called tristimulus values) are physically possible. Consider the LMS coordinates (0, 1, 0). To physically achieve this coordinate, we would need to find some way of stimulating the M cones without stimulating the L or S cones at all. Because the M cone’s sensitivity curve significantly overlaps at least one of L or S at all wavelengths, this is impossible!

A problematic side effect of this fact is that it’s really difficult to increase stimulation of only one of the cones. This, in particular, would make it not a great candidate for building display hardware.

Another historical, pragmatic problem was that the cone sensitivities weren’t accurately known until the 1990’s, and a need to develop a mathematically precise model of color significantly predates that. The first significant progress on that front came about in the late 1920’s.

Wright & Guild’s color matching experiments

In the late 1920’s, William David Wright and John Guild conducted experiments to precisely define color in terms of contributions from 3 specific wavelengths of light.

Even though they may not have known about the three classes of cones in the eye, the idea that all visible colors could be created as the combination of three colors had been proposed at least a hundred years earlier.

Wright & Guild had the idea to construct an apparatus that would allow test subjects to reconstruct a test color as the combination of three fixed wavelength light sources. The setup would’ve looked something like this:

The experimenter would set the lamp on the bottom to a target wavelength, (for instance, 600nm) then ask the test subject to adjust the three lamp power controls until the colors they were seeing matched.

The power settings of the three dials give us a (red, green, blue) triplet identifying the pure spectral color associated with 600nm. Repeating this process every 5nm with about 10 test subjects, a graph emerges showing the amounts of red (700nm), green (546nm), and blue (435nm) light needed to reconstruct the appearance of a given wavelength. These functions are known as color matching function (CMFs).

These particular color matching functions are known as $$\bar r(\lambda)$$, $$\bar g(\lambda)$$, and $$\bar b (\lambda)$$.

This gives the pure spectral color associated with 600nm an $$(R, G, B)$$ coordinate of (0.34, 0.062, 0.00). This is a value in the CIE 1931 RGB color space.

Hold on though — what does it mean when the functions go negative, like here?

The pure spectral color associated with 500nm has an $$(R, G, B)$$ coordinate of (-0.72, 0.85, 0.48). So what exactly does that -0.72 mean?

It turns out that no matter what you set the red (700nm) dial to, it will be impossible to match a light outputting at 500nm, regardless of the values of blue and green dials. You can, however, make the two sides match by adding red light to the bottom side.

The actual setup probably would’ve had a full set of 3 variable power, fixed wavelength lights on either side of the divider to allow any of them to be adjusted to go negative.

Using our color matching functions, we can match any monochromatic light using a combination of (possibly negative) amounts of red (700nm), green (546nm), and blue (435nm) light.

Just as we were able to use our L, M, and S cone sensitivity functions to determine cone excitation for any spectral distribution, we can do the same thing with our color matching functions. Let’s apply that to the lemon color from before:

By taking the area under the curve of the product of the spectral curve and the color matching functions, we’re left with an $$(R, G, B)$$ triplet (1.0, 0.8, 0.2) uniquely identifying this color.

While the $$(L, M, S)$$ color space gave us a precise way to identify colors, this $$(R, G, B)$$ color space gives us a precise way to reproduce colors. But, as we saw in the color matching functions, any colors with a negative $$(R, G, B)$$ coordinate can’t actually be reproduced.

But this graph only shows which spectral colors can’t be reproduced. What about non-spectral colors? Can I produce pink with an R, G, B combination? What about teal?

To answer these questions, we’ll need a better way of visualizing color space.

Visualizing color spaces & chromaticity

So far most of our graphs have put wavelength on the horizontal axis, and we’ve plotted multiple series to represent the other values of interest.

Instead, we could plot color as a function of $$(R, G, B)$$ or $$(L, M, S)$$. Let’s see what color plotted in 3D $$(R, G, B)$$ space looks like.

Cool! This gives us a visualization of a broader set of colors, not just the spectral colors of the rainbow.

A simple way to reduce this down to two dimensions would be to have a separate plot for each pair of values, like so:

In each of these plots, we discard one dimension by holding one thing constant. Rather than holding one of red, green, and blue constant, it would be really nice to have a plot showing all the colors of the rainbow & their combinations, while holding lightness constant.

Looking at the cube pictures again, we can see that (0, 0, 0) is black, and (1, 1, 1) is white.

What happens if we slice the cube diagonally across the plane containing $$(1, 0, 0)$$, $$(0, 1, 0)$$, and $$(0, 0, 1)$$?

This triangle slice of the cube has the property that $$R + G + B = 1$$, and we can use $$R + G + B$$ as a crude approximation of lightness. If we take a top-down view of this triangular slice, then we get this:

This two dimensional representation of color is called chromaticity. This particular kind is called rg chromaticity. Chromaticity gives us information about the ratio of the primary colors independent of the lightness.

This means we can have the same chromaticity at many different intensities.

We can even make a chromaticity graph where the intensity varies with r & g in order to maximize intensity while preserving the ratio between $$R$$, $$G$$, and $$B$$.

Chromaticity is a useful property of a color to consider because it stays constant as the intensity of a light source changes, so long as the light source retains the same spectral distribution. As you change the brightness of your screen, chromaticity is the thing that stays constant!

There are many different ways of dividing chromaticity into two dimensions. One of the common methods is used in both the HSL and HSV color spaces. Both color spaces split chromaticity into “hue” and “saturation”, like so:

It might appear at a glance that the rg chromaticity triangle and these hue vs. saturation squares contains every color of the rainbow. It’s time to revisit those pesky negative values in our color matching functions.

Gamuts and the spectral locus

If we take our color matching functions $$\bar r(\lambda)$$, $$\bar g(\lambda)$$, and $$\bar b(\lambda)$$ and use them to plot the rg chromaticities of the spectral colors, we end up with a plot like this:

The black curve with the colorful dots on it shows the chromaticities of all the pure spectral colors. The curve is called the spectral locus. The stars mark the wavelengths of the variable power test lamps used in the color matching experiments.

If we overlay our previous chromaticity triangles onto this chart, we’re left with this:

The area inside the spectral locus represents all of the chromaticities that are visible to humans. The checkerboard area represents chromaticities that humans can recognize, but that cannot be reproduced by any positive combination of 435nm, 546nm, and 700nm lights. From this diagram, we can see that we’re unable to reproduce any of the spectral colors between 435nm and 546 nm, which includes pure cyan.

The triangle on the right without the checkerboard is all of the chromaticities that can be reproduced by a positive combination. We call the area that can be reproduced the gamut of the color space.

Before we can finally return to hexcodes, we have one more color space we need to cover.

CIE XYZ color space

In 1931, the International Comission on Illumination convened and created two color spaces. The first was the RGB color space we’ve already discussed, which was created based on the results of Wright & Guild’s color matching experiments. The second was the XYZ color space.

One of the goals of the XYZ color space was to have positive values for all human visible colors, and therefore have all chromaticities fit in the range [0, 1] on both axes. To achieve this, a linear transformation of RGB space was carefully selected.

The analog of rg chromaticity for XYZ space is xy chromaticity and is the more standard coordinate system used for chromaticities diagrams.

Gamuts are typically represented by a triangle placed into an xy chromaticity diagram. For instance, here’s the gamut of CIE RGB again, this time in xy space.

With an understanding of gamuts & chromaticity, we can finally start to discuss how digital displays are able to display an intended color.

Screen subpixels

Regardless of the manufacturer of your display, if you took a powerful magnifying glass to your display, you would find a grid of pixels, where each pixel is composed of 3 types of subpixels: one type emitting red, one green, and one blue. It might look something like this:

Unlike the test lamps used in the color matching experiments, the subpixels do not emit monochromatic light. Each type of subpixel has its own spectral distribution, and these will vary from device to device.

Using ColorSync Utility on my Macbook Pro, I was able to determine the xy space gamut of my screen.

Notice that the corners of the gamut no longer lie along the spectral locus. This makes sense, since the subpixels do not emit pure monochromatic light. This gamut represents the full range of chromaticities that this monitor can faithfully reproduce.

While gamuts of monitors will vary, modern monitors should try to enclose a specific other gamut: sRGB.

sRGB

sRGB (“standard Red Green Blue”) is a color space created by HP and Microsoft in 1996 to help ensure that color data was being transferred faithfully between mediums.

The standard specifies the chromaticities of the red, green, and blue primaries.

If we plot these, we wind up with a gamut similar to, but slightly smaller than, the MacBook LCD screen.

There are parts of the official sRGB gamut that aren’t within the MacBook Pro LCD gamut, meaning that the LCD can’t faithfully reproduce them. To accommodate for that, my MacBook seems to use a modified sRGB gamut.

sRGB is the default color space used almost everywhere, and is the standard color space used by browsers (specified in the CSS standard). All of the diagrams in this blog post are in sRGB color space. That means that all colors outside of the sRGB gamut aren’t accurately reproduced in the diagrams in this post!

Which brings us, finally, to how colors are specified on the web.

sRGB hexcodes

#9B51E0 specifies a color in sRGB space. To convert it to its associated (R, G, B) coordinate, we divide each of the three components by 0xFF aka 255. In this case:

0x9B/0xFF = 0.61
0x51/0xFF = 0.32
0xE0/0xFF = 0.88

So the coordinate associated with #9BE1E0 is $$(0.61, 0.32, 0.88)$$.

Before we send these values to the display hardware to set subpixel intensities, there’s one more step: gamma correction.

Gamma correction

With each coordinate in RGB space being given 256 possible values, we want to ensure that each adjacent pair is as different as possible. For example, we want #030000 to be as different from #040000 as #F40000 is from #F50000.

Human vision is much more sensitive to small changes in low energy lights than small changes to high energy lights, so we want to allocate more of the 256 values to representing low energy values.

To see how, let’s imagine we wanted to encode greyscale values, and only had 3 bits to do it, giving us 8 possible values.

If we plot grey values as a linear function of energy, it would look something like this:

We’ll call our 3 bit encoded value $$Y$$. If our encoding scheme spaces out each value we encode evenly ($$Y = \frac{\left\lfloor8E\right\rfloor}{8}$$), then it would look like this:

You can see that the perceptual difference between $$Y=0$$ and $$Y=1$$ is significantly greater than the difference between $$Y=6$$ and $$Y=7$$.

Now let’s see what happens if we use a power function instead. Let’s try $$Y = \left(\frac{\left\lfloor8E\right\rfloor}{8}\right)^2$$.

We’re getting much closer to perceptual uniformity here, where each adjacent pair of values is as different as any other adjacent pair.

This process of taking energy values and mapping them to discrete values is called gamma encoding. The inverse operations (converting discrete values to energy values) is called gamma decoding.

In general form, gamma correction has the equation $$V_{out} = A V_{in}^\gamma$$. The exponent is the greek letter “gamma”, hence the name.

The encoding & decoding rules for sRGB use a similar idea, but slightly more complex.

If we plot sRGB values against linear values, it would look like this:

Okay! That was the last piece we needed to understand to see how we get from hex codes to eyeballs! Let’s do the walkthrough 😀

From Hexcodes to Eyeballs

First, we take #9B51E0, split it up into its R, G, B components, and normalize those components to be the range $$[0, 1]$$.

This gives us a coordinate of $$(0.61, 0.32, 0.88)$$ in sRGB space. Next, we take our sRGB components and convert them to linear values.

This gives us a coordinate $$(0.33, 0.08, 0.75)$$ in linear RGB space. These values are used to set the intensity of the subpixels on the screen.

The spectral distributions of the subpixels combine to a single spectral distribution for the whole pixel.

The electromagnetic radiation travels from the pixel through your cornea and hits your retina, exciting your 3 kinds of cones.

Putting it all together for a different color, we’re left with the image that opens this post!

A brief note about brightness setting

Before sRGB values are converted into subpixel brightness, they’ll be attenuated by the device’s brightness setting. So the 0xFF0000 on a display at 50% brightness might match the 0x7F0000 on the same display at 100% brightness.

In an ideal screen, this would mean that regardless of the brightness setting, black pixels $$(0, 0, 0)$$ would emit no light. Most phone & laptop screens are LCD screens, however, where each subpixel is a filter acting upon white light. This video is a great teardown of how LCDs work:

The filter is imperfect, so as brightness is increased, black pixels will emit light as the backlight bleeds through. OLED screens (like on the iPhone X and Pixel 2) don’t use a backlight, allowing them to have a consistent black independent of screen brightness.

Stuff I left out

This post intentionally glosses over many facets of color reproduction and recognition. For instance, we didn’t talk about what your brain does with the cone excitation information in the opponent-process theory or the effects of color constancy. We didn’t talk about additive color vs. subtractive color. We didn’t talk about color blindness. We didn’t talk about the difference between luminous flux, luminous intensity, luminance, illuminance, and luminous emittance. We didn’t talk about ICC device color profiles or what programs like f.lux do to color perception.

I left them out because this post is already way too long! As a friend of mine said: even if you’re a person who understands that most things are deeper than they look, color is way deeper than you would reasonably expect.

References

I spent an unusually large portion of the time writing this post just reading because I kept discovering that I was missing something I needed to explain as completely as I’d like.

Here’s a short list of the more helpful ones:

• A Beginner’s Guide to Colorimetry
• CIE 1931 Color Space
• HSL and HSV
• Gamma correction
• Real-Time Rendering, Third Edition p210-217

I also needed to draw upon many data tables to produce the charts in this post:

• University College London Color & Vision Research Laboratory database (XYZ color matching functions, cone fundamentals)
• fluxometer.com (RGB LCD screen subpixel spectra)
• CIE 15: Technical Report: Colorimetry, 3rd edition (RGB color matching funtions)

Special thanks to Chris Cooper and Ryan Kaplan for providing feedback on the draft of this post.

This is going to be a long post, so here’s the TL;DR that I want you to take away. I started experiencing my first bout of major depressive disorder in 2016, and am recovering now with the aide of meditation and 75mg daily of an anti-depressant called Zoloft (aka Sertraline). Talk therapy, medication, and meditation have all been important parts of my ongoing recovery. I’d love to talk to you if you think you might be in the same boat and want support pursuing help. If you don’t have the same symptoms I described here, please, please don’t take that to suggest you don’t need to seek some variety of help.

I recognize that those options are not available to everyone, and recognize that I’m fortunate to have the socioeconomic freedom to pursue all these options.

If you want something easier to digest, I’d suggest John Green’s post about his Mental Illness, and Hyperbole and a Half’s “Depression Part 2”. I hope this post can help a few people feel like they’re not alone in their struggle, just like those two posts did for me.

And because it will become blindingly apparent through this post, yes, indeed, I have gone full Valley: I’m an engineer in San Francisco working at a startup, attending meditation retreats, seeing mental health professionals, and went to Burning Man last year.

With that in mind, let’s go:

For most of my life, my emotional life has been pretty healthy. If I think about it in a ridiculously oversimplified good vs. bad feels over time framework, my emotional world might look like this:

At the beginning of 2016, my subjective world started to change.

The Onset

January - March 2016

The year started out on good footing, with motivation for work, exercise, living in a house of good friends, and finally having the chance to live with my girlfriend, Bonnie, again instead of mostly co-existing through walls of Facebook messenger chats and video calls across three timezones.

The first change seem innocuous enough. Standing outside the doorway of one of my housemate’s rooms, I commented that I’d noticed that I didn’t have the same visceral reaction to spending time with friends as I used to. He commented that it sounded pretty weird, and that maybe I should seek out therapy. I brushed off the notion, seeing it as a pretty minor problem, and went about my life. I started to notice the same thing would happen going on hikes. Any visceral sense of serenity was eerily missing.

My subjective world started to look like this:

The second change, though more gradual, was more obviously problematic. Over the course of months, I noticed my general mood decline. I’d have days waking up feeling kind of crummy with no obvious cause. As this evolved, my desire to interact with coworkers over lunch all but vanished. There were many days where I couldn’t bring myself to do it at all, so I hid in small single-chair rooms reserved for video calls.

The whole graph shifted down.

The last change was downright disturbing. I’d look down at my hands and experience them as just these objects I control. I’d look in the mirror and see a person that I rationally knew was me, but had no visceral sense of recognition. Photos of me in the past felt like looking at a bizarre alternative timeline that I could remember but not feel. From searching around, the best terms I’ve found to describe this are “dissociation” or “depersonalization”.

By the end of February, I knew something had to change.

So I did the first thing that enters any Valley yuppy’s head: quit my job to go travel the world to try to find myself. March 25th, 2016 was my last day at Khan Academy.

By March 31st, I had landed in Hong Kong on a one-way ticket.

Travel

April - June 2016

This is not a post about how travel helped me face my inner demons or expanded my comfort zone or created lifelong friendships. Travel has done wonders for me in the past, but it’s not exactly a silver bullet for mental health.

As my general mental health deteriorated further, a perverse reversal took place.

See that dip to the right of the green line? Yeah, that’s not normal, and it merits some explanation.

When I went to Hong Kong in 2014 with some friends, I remembered a distinct feeling of serenity standing in the Nan Liang Garden.

The sky, the trees, and the pagoda all brought me comfort, peace, and a gratitude for the world.

When I returned in 2016, here was my initial feeling: nothing. I might as well have been standing in the parking lot of Walmart.

So that reaction would explain a flatline, but why is there a dip? Let me show you an excerpt of the internal chatter of a chronic ruminator.

“Hey, I remember being here, last time it was super great.”

“Hmm, this time I don’t think I really feel the same connection.”

“Maybe I just need to relax more, let me give that a try.”

“Nope, still nothing. Well that’s disappointing.”

“Well if this doesn’t lift my mood, maybe nothing else will.”

“Why am I thinking this way, it’s self destructive — I should stop thinking this way.”

“I’ll bet this is going to happen every time I go anywhere nice that I’ve been before.”

“What’s wrong with me? Am I going to be like this forever?”

“What’s the point of living if I can only experience good things in a superficial way?”

Positive experiences had now become an exercise in “let’s painstakingly analyze how deeply broken my emotional mind is”. Not the most constructive internal dialogue.

Every time anything good would happen, I would cling on as hard as I could to any positive aspect of it, then be intensely distraught when it would inevitably fade almost immediately after I had left the experience.

Perhaps the most soul-crushing parts of the trip were the periods where I actually did feel better for days at a time. Each time it happened, I had this gloriously optimistic thought: “Maybe it’s over. Maybe I’m going to feel normal again from now on!” And each time I could feel my sensitivity to the world receding, it felt like a force pulling away everything good in my present and future. I could never figure out a way to return to that feeling of normalcy. Once, I felt that lucidity kick in for about 10 minutes during a car ride, then mysteriously vanish right afterwards.

Travel had, in previous long trips, been about novel experiences and collecting stories. But now novel experiences were just one more opportunity to scrutinize my reduced sensitivity. And collecting stories isn’t much fun when they don’t feel like your own stories.

So, after hopping on and off planes for a few months, I returned to California to look for work and seek professional help.

Therapy & Tears

July - September 2016

I’m fortunate to have had enough friends, family, and acquaintances share their experiences going to see mental health professionals that I had a significantly reduced mental block for pursuing help. I hope I can take one more brick off that wall for those reading this.

On advice from my friend David, I emailed four different therapists asking for a free 15 minute phone consultation (these consultations are available from most therapists). They replied to set up times, and I found that for the half hour leading up to each calls, I was basically non-functional and reverted to lying in bed watching Key & Peele videos on YouTube to curb my nervousness.

During the first of the calls, I described my problems, and the person very politely said that yes, they thought it was something they could probably help me with. That seemed like a positive signal to me, but I figured I’d wait to hear from the other therapists.

On my next call, I spoke to a therapist near Redwood City named Sandy Lillie. I described my problems, and she very matter-of-factly told me it sounded like I had been through something traumatic and had emotionally shut down.

Unless you’ve been there yourself, it’s hard to appreciate how validating it is to hear from an authoritative source that yes, you do have a real problem. It’s not just “part of getting older” or a “slump”, it’s a real problem that can and should be addressed.

That validation was so powerful that I spent the next half hour as a teary mass under the blankets in my room. I cancelled the rest of the phone consultations and booked my first 75 minute appointment.

During the first appointment, I told Sandy about my relationship with Bonnie, my family, and a bit about how I view myself and work.

I cried. A lot. Like half a waste basket full of tissues a lot. The most surprising part to me was not that I cried, but when I cried. Areas of my life I had never considered as harbors of pain were being poked and prodded and all their juicy bits were spilling out of my eyeballs into the ample supply of Kleenex provided.

I’m grateful at this point that I’ve been keeping journals for years and occasionally taking the time to put a stream-of-consciousness dump in there, because I find it remarkably difficult to remember how I felt at a given time. Thankfully, I can have July 2016-me tell the story of driving home after my first appointment:

When I got in the car and started driving back, I found myself in this bizarre state between laughing and crying half of the way back. Then I started worrying that we had a flat and was running low on gas so my brain started focusing on practicalities again. For a little while there, it almost felt like the symptoms of dissociation were fading, and that I could experience the beauty of the landscape I was driving through.

I’m going back again on Monday. That feeling of laughing and crying is incredible.

After every few sessions of therapy, I’d have a huge emotional release. It felt like there were these buried pressurized tanks of feelings that were being dug up and punctured. Whenever it happened, I felt genuinely better in a sustained way for a few days.

When I noticed this trend, I started to go hunt for things in my past for emotional release outside of therapy sessions. I read journal entries from up to 4 years past to try to extract some nourishing tears. Perversely, whenever I noticed myself feeling especially dead, I’d go hunting for ways to make myself cry. It really was true that feeling something other than “dead” or self-pitying was significant.

Every time I did this though, the relief period got shorter and less intense.

Until I ran out of past sadness to extract.

The only thing I would cry about was how dead I felt, which, unlike digging up my past, did not feel cathartic at all. It felt pathetic. Useless. Shameful.

Feeling Supported

August 2016

At that point, I started feeling deeply disconnected again. I realized that my visceral connection with Bonnie was all but gone, and that she could feel it too. I told her that when I looked into her eyes, I wanted desperately to feel connected.

People talk about eyes acting as the windows into people’s soul. Through them, you see your past experiences with them, your future plans, your shared hopes, and the struggles you’ve overcome together.

But all I saw was glass.

It was the difference between reading a love letter and staring at illegible pen scratches on a piece of dead tree. The difference between your home and a house staged for sale. A baseball caught from a home run vs one bought at the store. In each case, each item of the pair are objectively the same, but what they each represent is fundamentally different.

We talked about our future, and she shared her vision for it. I admitted I couldn’t see it. She asked if I wanted to keep trying. I said yes. She asked why. I struggled.

Here’s August 2016-me’s account of what followed:

I squirmed uncomfortably, feeling dread that I didn’t have an answer, until it rushed back with “I want to go back to Kitchener”. And I sobbed. Hard. [Kitchener was the first place we’d lived together — we shared an apartment for 11 months starting in 2014.]

After I said that, and kept crying, I felt like I could connect to my past. I felt like me. I felt like I found myself underneath this judgmental voice in my head. I talked about good times in Kitchener. Cooking together, going to the grocery story together, rock climbing together. I told her I just wanted the simple life back, and wanted to be a kid.

…

I could feel it go in and out. I could feel the true part of me come back to me, then recede behind whatever’s blocking it. When it was there, memories of the past flooded back to me, and they felt like mine. I was able to recognize myself in the past. I told Bonnie that everything felt smaller when I felt like myself. Like we were kids.

I told her I was afraid we couldn’t go back. But that I wanted to try. And I was afraid that the morning after, I wouldn’t be able to get back to this. Bonnie told me that she thinks she does get a visceral reaction from looking into my eyes, but hasn’t seen anything there in a while.

After using up nearly all the toilet paper for our tears, and after I felt like the real me had receded back behind the barrier entirely for the night, we crawled into bed and went back to sleep.

The next day I felt dead again, but from then on I knew that Bonnie understood there really was something wrong with me. I had support.

Someone close to me finally knew I was honestly trying.

Therapy & Pure Rage

September 2016

Returning from the emotional overwhelm that is Burning Man, I had a new insight into my emotional life: I am absolutely awful at productively expressing anger.

Sometime during high school, it became apparent to me that my default reaction mode while angry was to lash out, and that was less than ideal. As the introspective voice in my head expanded its domain over my life, at some point anger came under its purview. The voice’s approach was to try to explain to me why anger was not the productive response to situations, and why my anger was irrational and unfair most of the time.

To some extent, this is a useful skill, but I didn’t realize just how much I had bottled up until my fifth therapy session.

During that session, Sandy asked if I could role play talking to a specific person that had deeply hurt me. I agreed, and as I started to express myself coldly, Sandy pushed me to be more expressive.

Here’s September 2016-me recounting what followed:

I could feel my whole body trembling with rage, and could feel myself fighting internally to express it. Sandy tried to encourage me to put some volume behind “Fuck You”, but I found it was too hard for me. Instead, I said it quietly between tears. Eventually, I managed to hit a threshold, and started screaming. No words, just raw, vocally expressed anger. I think I did that for about half an hour, almost until the end of the session. When I opened my eyes, the world felt a little bit more vivid again, like it did after the first time I had a long crying session during therapy.

It’s hard to describe the sound I made during that session. I had definitely never heard that sound come out of a human, even in movies. It didn’t feel so much like something I was doing as it did like a compressed demonic gas rapidly exiting my throat, dragging claws along my vocal chords on the way out. It wasn’t the last time I heard that sound.

The experience of the following 2 weeks or so was predictable:

By November, pretty much nothing seemed to be helping, and I’d begun to experience significant social anxiety even when I was in groups exclusively composed of friends I’d known for years. My general mood decline had also steadily continued, with probably a day or two a week starting out in bed crying for no discernible reason.

When I say “for no discernible reason”, I don’t mean that I had some minor bad experience that I was blowing out of proportion, I mean that between sobs I was contemplating what was externally wrong and couldn’t think of a single thing.

Have you ever tried to comfort someone who feels like crap for no reason? It’s not such an easy task.

Meditation

November 2016 - Present

For American Thanksgiving, I had a friend, also coincidentally named Sandy, over for dinner. Sandy was the second of my friends to become a dedicated meditator, and had spent a lot of time digesting hours and hours and hours of meditation audiobooks. He’d seen such dramatic improvements in his life from meditation that he was getting interested in trying to share some of those benefits with friends (or honestly, with whoever would listen).

When I told him how I was feeling and he suggested meditation, I was initially a little hesitant. I’d meditated fairly regularly using the Headspace app at the end of 2015, but had stopped once depression had set in for two reasons. The first was that I felt even more disconnected from the world after each meditation session than before, and the second was that when my mind felt scattered, the meditation sessions just seemed to repeatedly emphasize to me how little control over my brain I did have.

Sandy helped me a lot with the second objection by reframing meditation for me. He told me that the immediate goal during meditation is not long-held focus, it’s to notice when it drifts away. This sounds like a trivial difference, but it was key in how I responded to my mind drifting.

Each time my mind drifted before, I’d feel like I was failing and would be upset by how scattered my brain was. After the reframe, each time I’d notice it drift, I’d have this feeling of “Aha! I noticed! I succeeded!”.

With this framing, the world became just slightly more vivid after each session. It was still frustrating a significant amount of the time both during and after the session, but it paid off in the end. My mood began to stabilize back to where it was a few months earlier.

After meditating, I still wouldn’t generally feel like seeking out social interaction, but at least I didn’t dread it.

In December 2016, I had what I’d consider my first major mood reversal due to meditation. The context was that I was feeling terrible and having a hard time finding a quiet place to meditate because of roommates. Eventually, I locked myself in the bathroom with the fan on to drown out any sound and sat down to meditate on the toilet.

Here’s my journaled account of what happened that night during meditation:

Get very frustrated at night and break down and cry, but manage to break the cycle! Another part of my mind broke the downward spiral by repeating “You didn’t give up!” over and over and over again until the tears stopped and I felt happy again. I looked at face covered in tears and snot in the mirror, and was okay with what I saw :)

Despite the occasional boost from meditation, I was still feeling generally horrible on the whole. When I went back home for Christmas, I discovered that my social anxiety extended not only to friends, but also to family. I found myself frequently wanting to exit conversation with my parents and barricade myself in my room. I spent a lot of time using my terrible coping mechanism of burying myself in information through page after page of reading Reddit.

Anti-Depressants

January 2017 - Present

By January 2017, I turned to the option I told myself I would only turn towards if nothing else had worked: medication.

After the last visit to my therapist, she told me that from all the things she had heard from me, she could understand why I felt pain, but not of this magnitude. The struggles I had were real, but not uncommon. She suggested that it might be wise to try medication, and told me I’d need to go see an MD for that since she couldn’t prescribe medication (she has a PhD in Clinical Psychology).

From contemplating it for a year and going through therapy, I was able to precisely articulate what my emotional world looked like, and in two sessions with a psychiatrist, I was diagnosed with Major Depressive Disorder. He offered many potential treatment paths, among them medication.

On January 20, 2017, he prescribed me 50mg daily of Zoloft, with a 25mg daily ramp up period for a week. He told me that I shouldn’t expect to see any benefit in the first two weeks, and that was just a trial period meant to evaluate whether the side effects were tolerable.

Despite seeing medication as a reasonable option for other people, the first day of starting medication still felt incredibly hard. In part, it was because I was afraid of side effects like suicide ideation. Thankfully, suicidal ideation and urges have not been part of my particular flavor of depression.

The other part of the internal struggle was that it felt, in some ways, like admitting that I was entirely broken, and had failed to heal myself because I was weak. “Weak and useless” were a recurring pair of words in my journals and internal dialogue for most of 2016.

I don’t think I would’ve had the courage to start if it were not for support from family and for reading John Green’s post about his mental illness.

I know almost the exact point where the medication started working. It was late at night on January 26, 2017.

I was feeling pretty shitty through dance and afterwards, but decided I should meditate at home. After I did, I felt better. Near the end of the meditation, I wanted to give up, but didn’t let myself. Just as I previously had a voice tell me “You didn’t give up”, this time I had a voice tell me “You will get stronger” and repeat that over and over again.

The day after that, I felt better too. And the day after that.

I’d been burned by previous desperate optimism, so I was hesitant to believe this was truly the beginning to a sustained period of emotional stability.

But now, four months later, I can fairly confidently proclaim: I’m stable.

Interestingly, one of the first things to vanish as the result of the medication taking effect was any aversion to taking the medicine. Seeing it as a personal failing that I couldn’t work through myself started seeming as ridiculous as having an aversion to eye glasses because it demonstrates a failure on your part to will your eyeballs to reshape themselves.

In 2013, there were more than 41 million prescriptions of Zoloft alone, and that’s just in the US! That’s more than 1 for every 10 people in the country! What a bizarre thing to feel a social taboo over.

That said, there is something unsettling about the fact that after a year of struggle, tears, disconnection, existential crisis, and generally shitty feelings, the catalyst for recovery is half the size of a penny and costs me about $25/mo after insurance.

Brain chemistry is a hell of a thing.

Restarting Life

February 2017

Once I became more confident I was emerging from depression, I realized that I had to actually make some choices in life again. Up until that point, my decisions about what to do were largely guided by “what is the least likely to cause catastrophic mental meltdown?”.

Trying to change any of my negative habits was totally out of the question, since failing to do that would’ve been a crushing sign of additional failure in a life that felt full of weakness and uselessness.

But once I started recovering, I had a newfound certainty that for most plausible things that could happen in the external world, my internal world would be okay. So I started to make changes.

The first thing to go was the emotional crutch I’ve held for years of drowning myself in information to plaster over things I don’t want to think about. I quit my Facebook news feed, Twitter, Hacker News, and Reddit cold turkey. More generally the rule is to stay off sites that are infinite time sinks (which usually have “engagement” as their core metric). “How Technology is Hijacking Your Mind” is a great overview of how horribly misaligned many technology products are with the best interests of their users.

With the exception of a few days checking Twitter while at a conference, I’ve stayed completely off since February 5th. The time was replaced by watching educational content on YouTube, and reading books. When I get stuck at work, I just get up and walk around, maybe grab a healthy snack, maybe chat with coworkers.

It was interesting to see that in the first three weeks or so, I still had a compulsive habit of unlocking my phone and searching for those applications. I would opening Chrome and start typing “reddit.com” or “news.ycombinator.com”, before catching myself and closing it. I did this several times a day, which demonstrates oh-so-clearly to me how much of a compulsion it was and how little rational decision making went into the use of those websites.

Even now, four months later, I occasionally find myself sitting in front of my computer mentally cycling through all the sites I used to frequent to figure out which one I haven’t looked at in the last hour. Thankfully, I usually now have the sense to just close my laptop when I notice this happening. I suspect I’m much more aware of when this is happening due to the mindfulness meditation practices I’ve been doing.

The second thing was to organize my goals revolving around a monthly theme. Last month was “friends” — I enjoyed spending time with friends again! I decided to invite friends over to cook with, and reached out to some friends to video chat, one of whom I hadn’t talked to in years. It was great.

Meditation May

May 2017

May’s theme, now just finishing up, was meditation.

You might be wondering why I was still pursuing meditation when it had seemingly fallen far short of the improvement seen via pharmaceuticals.

Well, for starters, I consider the pills to be a catalyst for stability, but I suspect meditation had a lot to do with the rate of improvement. I started seeing benefits in the first two weeks, when I was told to not expect much until week 4-6.

Secondly, I don’t want to take pills forever. The recommended period is to stay with the dosage for a year after the depressive symptoms subside, and I’d like to do whatever I can to prevent relapse once I stop medicating. That said, if I had to choose medicating forever vs. dealing with depressive symptoms forever, I’d choose medication hands down. I’d just rather do what I can to deal with neither ☺️. Mindfulness meditation seems to be one of the few things with statistically significant results in preventing remission. One study found that it was as effective as maintenance level dosages of anti-depressants for preventing remission.

Lastly, and perhaps most critically, the pills boosted and stabilized my baseline mood, but things still felt a bit off.

Before the pills, I would feel shitty for no apparent reason. Afterwards, I would feel pretty-okay for no apparent reason. Pleasant experiences no longer had the perverse negative effect they previously had, but they also didn’t deliver much by way of connection or lasting joy.

But why would I expect meditation to help with this? Isn’t meditation only helpful for de-stressing or halting incessant mental chatter?

To explain, let’s talk a bit about lunchtime.

Eating Mindfully

Let me describe the process of eating as it routinely took place for me in years past.

Once seated to eat, I’d collect a spoonful of food, and, if I’m lucky, stay focused long enough to taste the first bit going into my mouth. For each bite after that, some part of my brain would be scheduling my day or planning for the future or replaying an awkward encounter from the previous day. Another part of my brain would be guiding my eyes to provide enough visual information for my hand to construct a continuous stream of sustenance to be delivered to my mouth, as my teeth and throat work in tandem to get those nutrients to my stomach. That leaves just about…. zero percent of my conscious brain to actually taste the damn food. Unsurprisingly, this meant that most meals felt unremarkable.

Now I try to do the following: collect a spoonful of food, look at it carefully, focus visually on what I’m about to eat. Next, put it in my mouth, and focus on the taste and texture of the food, ideally with my eyes closed. Once I’ve fully tasted it and swallowed, only then do I go in search of my next bit of food to eat from my plate and return to step one.

This might seem silly or trivial, but give it a try on your next meal, even if it’s just your morning toast. You might discover how little you’ve been tasting your food too.

So what does this have to do with meditation? Everything! Meditation can systematically cultivate your capacity to concentrate on the most presently relevant thing. The critical bit needed to enjoy the taste of your food is the ability to repeatedly and gently return your attention to the thing that’s most important in the current moment.

In the book “The Mindful Way Through Depression”, the authors describe the experience of one depressed individual mindfully eating a single raisin:

“I was very aware of what we were doing. I’ve never tasted a raisin like that before. Actually, I’ve never noticed what a raisin looked like. At first it looked dead and crinkly, but then I noticed how the light struck it in different ways, like a jewel. When I first put it in my mouth, it was hard at first to stop myself from instantly chewing it. Then, when I was exploring it with my tongue, I was able to tell which side was which — but there was no taste. Then when eventually I bit down on it — wow, it was absolutely amazing. I’d never tasted anything like it.”

During his lectures at the meditation retreat I recently returned from, our teacher, Shinzen Young, recounted the story of how he first knew his meditation was bearing fruit. Every morning at the monastery in Japan, their breakfast would be the same thing: plain rice porridge with no seasoning. One day he woke up and realized he was excited for his morning rice porridge. Shinzen grew up in California, and yet he couldn’t remember being as excited for belgian waffles with syrup, fresh hash browns, and fresh squeezed orange juice as he was for this bowl of rice and water.

Through meditation, it’s possible to systematically cultivate control over your attention in ways that will allow you appreciate everything in life in a deeper way.

Return of Sensitivity

After meditation, I notice that everything just seems more interesting. The position of birds wings as they swooped, the complex cacophony of sounds everywhere I go, the feeling of my bedsheets on my skin, the feeling of my shoe when I first slip it on, even the stucco pattern of the walls in my apartment.

As I’ve racked up hours of meditation over the last few months, I’ve noticed that the interesting-ness of everything tends to be present to a small degree even when I haven’t meditated in a while. It’s always subtly there in the background.

An even stranger phenomenon that I haven’t quite reached, but that I have friends that can attest to is this:

I won’t say too much about this, but suffice to say I have a friend that now actively seeks out opportunities to eat uncomfortably hot peppers by themselves, or to sit for a while without a shirt in a snowstorm, and emerge from it quite happy. Personally, I’ve come to see that I can quite enjoy freezing cold showers.

What Didn’t Work for Me

I tried a number of things that I didn’t dedicate sections in this post to because they didn’t do much good for me.

Exercise

Exercise, while I’m sure good for my muscles, didn’t have a clear positive impact on my mental health. On days where I didn’t perform as well at the gym as I’d hoped I would, it generally felt net-negative. The evidence for exercise combatting depression does seem pretty strong even when compared to Zoloft/Sertraline specifically, but that may be because the study group was starting from a baseline of relative inactivity, whereas I think I was already comfortably getting more than the recommended 30 minutes of exercise 3 times weekly.

The act of scheduling exercise on a regular basis and setting goals was, however, quite helpful. I would plan out my entire day hour by hour to prevent decision fatigue which could bring my day to a lying-on-the-floor-staring-at-the-ceiling halt.

Here’s one such schedule from my journal:

7:30am wakeup, brush teeth, smoothie
8:00am drive to gym
8:30am gym
9:45am drive home
10:00am shower
10:10am 10min headspace
10:20am setup time machine hard drive downstairs
10:30am watch ML while eating peanut butter toast
12:00pm lunch, eating lentil soup, continue watching ML
12:30pm 10min headspace
12:40pm Work on signed distance fields blogpost
5:00pm Dinner: lazy fried rice
6:00pm Walk to MV caltrain

Supplements

I tried taking high EPA fish oil and vitamin D supplements since they’ve both been somewhat effective for treating depression in studies. I didn’t notice any significant effect, and the fish oil was especially uncomfortable for me since I’m vegetarian, so I stopped. Some of the studies suggested that the fish oil was only effective for people already on anti-depressants, but I had already stopped taking the fish oil once I started on Zoloft.

Talking to Friends

I really wish that I found talking to friends about depression helpful, but I can’t honestly say it was. This was really only helpful to the extent that I felt I could help others seek help for their own conditions and to the extent that it encouraged me to seek help for myself.

There’s a special kind of frustration in knowing that you have friends that want to help you and not knowing a single thing they can do to help. I found that talking about my depression generally left me feeling self-pity afterwards, and not at all more understood, despite the earnest efforts of many friends.

Learned Optimism

I sought out the book “Learned Optimism” after it was described in a conference talk by Reginald Braithwaite, but quickly found that the description of depression in that book didn’t align with what I was experiencing. I found the descriptions in “The Mindful Way through Depression” much more relatable.

Closing Thoughts

So that’s my journey so far. Come next year as I try to step away from anti-depressants, I might have more to say. I consider myself quite fortunate to have been so responsive to pharmaceutical treatments at relatively low dosages with no noticeable side effects, and recognize that it’s not smooth sailing for everyone that chooses to pursue that path.

If you’re interested in reading more about meditation, I recommend Shinzen Young’s “See, Hear, Feel” introduction. It’s long, but pretty information dense. If you can’t be bothered to invest so much time up front, Headspace has a low initial commitment. If you’re interested in retreats, you can find the schedule of his retreats here: Retreats - Shinzen Young. I appreciate that he doesn’t ask you to dogmatically believe much of anything — he’s largely decoupled meditation practices from their religious doctrines, and hasn’t sacrificed depth of knowledge in the process.

If you want to talk to me about my experiences with medication, therapy, mindfulness, or share your own, you can reach by email at jamie.lf.wong@gmail.com (or Facebook messenger if we’re friends).

Best of luck out there in the world — it’s a tricky place!

This post is also available in Russian: Симуляция физического мира.

How might you go about simulating rain? Or any physical process over time, for that matter?

In a simulation, whether it be rain, the airflow over a plane wing, or a slinky slinking down some stairs, we can frame the entire simulation with two pieces of information:

1. What is the state of everything at the beginning of the simulation?
2. How does that state change from one moment of time to the next?

By “state of everything”, I mean any non-constant information needed either to determine what the scene looks like at a given moment, or about how the scene will change from one moment to the next. The position of a raindrop, the direction of the wind, and the velocity of each piece of the slinky are all examples of state.

If we represent the entire state of our scene with one big vector \( \vec y \), then we can mathematically reformulate the two pieces of information above into these two statements:

1. What is a value \(y_0\) that satisfies \( y(t_0) = y_0 \)?
2. What is a function \(f \) that satisfies \( \frac{dy(t)}{dt} = f(t, y(t)) \)?

If you’re confused about why you’d want to store the whole state in one big vector, bear with me. This is one of those cases where it might seem we go way over the top with generality, but I promise there are some interesting tidbits from looking at a general model for simulation.

If you’re confused about how you might store the whole state of a scene in one big vector, let’s look at a simple example. Let’s consider a 2D simulation with 2 particles. Each particle has a position \( \vec x \) and a velocity \( \vec v \).

So to make \( \vec y \), all we have to do is smoosh \( \vec x_1 \), \( \vec v_1 \), \( \vec x_2 \), and \( \vec v_2 \) into one 8 element vector, like this:

And if you’re confused why we want to find \( f(t, y(t)) \) instead of any old definition of \( \frac{dy(t)}{dt} \), the point is to find the derivative only terms of our current state, \( y(t) \), constants, and the time itself. If that’s impossible, it likely that there’s some bit of state we forgot to consider.

The Initial State

Defining \( \vec y_0 \) defines the initial state of the simulation. So if the initial state of our two-particle simulation looks something like this:

Then our numerical initial conditions might look like this:

Smooshing that together into a single vector, we get our \( \vec y_0 \).

The Derivative Function

\( \vec y_0 \) tells us the initial state, so now all we need is some way of getting from the initial state to the state a tiny bit into the future, and from that state a tiny bit further, and so on.

With that in mind, let’s solve for \( f \) in the equation \( \frac{dy(t)}{dt} = f(t, y(t)) \). So let’s take the derivative of \( y(t) \).

Woah! That’s a tall formula! Don’t worry, we can hopefully make it feel less intimidating by breaking \( \vec y \) back out to its constituent parts.

\( \vec x_1 \) and \( \vec x_2 \) are going to follow similar rules to one another, as will \( \vec v_1 \) and \( \vec v_2 \). So despite the ball of notation above, all we really want to find are the following two things:

How we define those two derivatives is the real nuts and bolts of the simulation, and to make this an interesting simulation and not just a program where random stuff happens, we can to look to physics for inspiration.

Kinematics & Dynamics

A bit of introductory kinematics and dynamics will go a long way in making our simulation interesting. Let’s start with the real basics.

\( \vec x \) represents position, and the first derivative of position with respect to time is velocity, \( \vec v \). In turn, the first derivative of velocity with respect to time is acceleration, \( \vec a \).

It might seem by now that we’ve already answered our question of finding our derivative function \( f \), since we know the following:

We have indeed rather nailed down \( \frac{d \vec x}{dt} \) since \( \vec v \) is part of our state vector \( \vec y(t) \), but we’ll need to go a teensy bit further for the second equation since \( \vec a \) is not.

Newton’s second law comes in handy here: \( \vec F = m \vec a \). If we assume the mass of our particles are known, then we can re-arrange that equation and we’re left with this:

Well hold on, \( \vec a \) wasn’t part of \( \vec y(t) \), neither is \( \vec F \), so this hardly seems like progress (remember, we need our derivative function to only be a function of \( \vec y(t) \) and \( t \)). But indeed it is progress, because we have all sorts of handy formulas that govern forces in the natural world.

Let’s pretend, for our simple example, that the only force acting upon the particles is the gravitational attraction between them. In that case, we can determine \( \vec F \) using Newton’s law of universal gravitation:

Where \( G \) is the gravitational constant \( 6.67 \times 10^{-11} \frac{Nm^2}{kg^2} \), and \( m_1 \) and \( m_2 \) are the masses of our particles (which we assume are constant).

In order to do our simulation, we need a direction too, and we also need to define \( r \) in terms of some part of \( \vec y(t) \). If we say that \( \vec F_1 \) is the force acting upon particle 1, then we can do that like so:

Let’s recap. The changing state of our two particle system is entirely described by \( \vec x_1, \) \( \vec v_1, \) \( \vec x_2, \) and \( \vec v_2 \). And those change over time like so:

Now we have all the information that makes this simulation different from all other simulations: \(\vec y_0\) and \(f\).

Now then, now that we have a rigorously defined simulation, how might we go about turning it into a mesmerizing animation?

In case you’ve written simulations or games before, you might jump to something like this:

x += v * delta_t
v += F/m * delta_t

Let’s take a step back and consider why that works.

A Differential Equation

Before diving into implementation, let’s step back a second and see what information we have and what information we need. We have a \( y_0 \) that satisfies \( y(t_0) = y_0 \), and we have an \( f \) that satisfies \( \frac{dy}{dt}(t) = f(t, y(t)) \). What we want is a function that can predict the state of the system at any point in time. Phrased mathematically, we want \( y(t) \).

With this in mind, if you squint carefully at \( \frac{dy}{dt}(t) = f(t, y(t)) \), you might spy that it’s an equation that relates \( y \) with its derivative \( \frac{dy}{dt} \). That makes it a differential equation! More specifically, it makes it a first order, ordinary differential equation. If we solve this differential equation, we get the function \( y(t) \).

Solving for \( y(t) \) given \( y_0 \) and \( f \) is called an initial value problem.

Numerical Integration

Some initial value problems have easily found analytic solutions, but for complex simulations this tends to be infeasible. So let’s try to figure out a method of approximating the solution.

Let’s explore a simpler initial value problem:

Given \( y(0) = 1 \) and \( \frac{dy}{dt}(t) = f(t, y(t)) = y(t) \), find an approximation of \( y(t) \).

Let’s examine the problem from a graphical perspective by considering the value and tangent line at \(t = 0 \). From our givens, we have \( y(0) = 1 \) and \( \frac{dy}{dt}(t) = y(t) = 1 \).

So we don’t know what \( y(t) \) looks like yet, but we do know that it should follow the tangent line near to \( t = 0 \). So let’s estimate \( y(0 + h) \) for some small value of \( h \) by following the tangent line. We’ll use \( h = 0.5 \) for now.

Symbolically, we’re estimating \( y(h) \) like so:

So for \( h = 0.5 \), \( y(h) \approx 1.5 \).

Now we can repeat this process. Even though we don’t know the exact value of \( y(h) \), as long as we have a pretty good approximation of it, we can figure out a pretty good approximation of the tangent line at that point too!

Then we can follow this tangent line \( h \) units to the right as well.

We can repeat the process again, with a tangent slope of \( f(t, y(t)) = f(1, 2.25) = 2.25\):

Stating the procedure recursively, we have this:

This is called the “Forward Euler” method of numerical integration. This is the general case of the step x += v * delta_t!

In this particular initial value problem, our steps look like so:

This method lends itself well to representing the computations in a table, like so:

It turns out that this particular initial value problem has a nice analytical solution: \( y(t) = e^t \)

When approximating the solution with Forward Euler, what do you think happens as the time step \( h \) gets smaller?

The error between the approximate solution and the exact solution decreases as \( h \) decreases! In addition to decreasing \( h \) to decrease the error, you could also use an alternative method of numerical integration that may provide better error bounds, such as the Midpoint method, Runge-Kutta methods, and linear multistep methods.

Let’s get coding!

Just as we can generalize the definition of a simulation mathematically, we can generalize the implementation of a simulation programmatically.

Since I’m most familiar with JavaScript but appreciate the clarity type annotations can give code, all the code samples will be written in TypeScript.

Let’s start with a version that assumes that y is a one-dimensional array of numbers, just as in the mathematical explorations.

function runSimulation(
    // y(0) = y0
    y0: number[],
    // dy/dt(t) = f(t, y(t))
    f: (t: number, y: number[]) => number[],
    // display the current state of the simulation
    render: (y: number[]) => void
) {
    // Step forward 1/60th of a second at a time
    // This results in realtime simulation if playback is 60fps
    const h = 1 / 60.0;

    function simulationStep(ti: number, yi: T) {
        render(yi)
        requestAnimationFrame(function() {
            const fi = f(ti, yi)

            // t_{i+1} = t_i + h
            const tNext = ti + h

            // y_{i+1} = y_i + h f(t_i, y_i)
            const yNext = []
            for (let j = 0; j < y.length; j++) {
                yNext.push(yi[j] + h * fi[j]);
            }

            simulationStep(tNext, yNext)
        })
    }
    simulationStep(0, y0)
}

It turns out to be pretty inconvenient to always operate on 1-dimensional arrays of numbers, so we can abstract away the addition and multiplication operations on the simulation state into an interface, then define our general simulation code concisely using TypeScript Generics¹.

interface Numeric<T> {
    plus(other: T): T
    times(scalar: number): T
}

function runSimulation<T extends Numeric<T>>(
  y0: T,
  f: (t: number, y: T) => T,
  render: (y: T) => void
) {
    const h = 1 / 60.0;

    function simulationStep(ti: number, yi: T) {
        render(yi)
        requestAnimationFrame(function() {
            // t_{i+1} = t_i + h
            const tNext = ti + h
            // y_{i+1} = y_i + h f(t_i, y_i)
            const yNext = yi.plus(f(ti, yi).times(h))
            simulationStep(yNext, tNext)
        })
    }
    simulationStep(y0, 0.0)
}

The cool thing about having this general code is it lets us focus on the core of the simulation: what about this simulation is different from other simulations. Using the example of our two particle simulation from before:

// Represents the state of the entire
// two particle simulation at a point in time.
class TwoParticles implements Numeric<TwoParticles> {
    constructor(
        readonly x1: Vec2, readonly v1: Vec2,
        readonly x2: Vec2, readonly v2: Vec2
    ) { }

    plus(other: TwoParticles) {
        return new TwoParticles(
            this.x1.plus(other.x1), this.v1.plus(other.v1),
            this.x2.plus(other.x2), this.v2.plus(other.v2)
        );
    }

    times(scalar: number) {
        return new TwoParticles(
            this.x1.times(scalar), this.v1.times(scalar),
            this.x2.times(scalar), this.v2.times(scalar)
        )
    }
}

// dy/dt (t) = f(t, y(t))
function f(t: number, y: TwoParticles) {
    const { x1, v1, x2, v2 } = y;
    return new TwoParticles(
        // dx1/dt = v1
        v1,
        // dv1/dt = G*m2*(x2-x1)/|x2-x1|^3
        x2.minus(x1).times(G * m2 / Math.pow(x2.minus(x1).length(), 3)),
        // dx2/dt = v2
        v2,
        // dv2/dt = G*m1*(x1-x1)/|x1-x2|^3
        x1.minus(x2).times(G * m1 / Math.pow(x1.minus(x2).length(), 3))
    )
}

// y(0) = y0
const y0 = new TwoParticles(
    /* x1 */ new Vec2(2, 3),
    /* v1 */ new Vec2(1, 0),
    /* x2 */ new Vec2(4, 1),
    /* v2 */ new Vec2(-1, 0)
)

const canvas = document.createElement("canvas")
canvas.width = 400;
canvas.height = 400;
const ctx = canvas.getContext("2d")!;
document.body.appendChild(canvas);

// Display the current state of the simulation
function render(y: TwoParticles) {
    const { x1, x2 } = y;
    ctx.fillStyle = "white";
    ctx.fillRect(0, 0, 400, 400);

    ctx.fillStyle = "black";
    ctx.beginPath();
    ctx.ellipse(x1.x*50 + 200, x1.y*50 + 200, 15, 15, 0, 0, 2 * Math.PI);
    ctx.fill();

    ctx.fillStyle = "red";
    ctx.beginPath();
    ctx.ellipse(x2.x*50 + 200, x2.y*50 + 200, 30, 30, 0, 0, 2 * Math.PI);
    ctx.fill();
}

// Run the thing!
runSimulation(y0, f, render)

If we tweak the numbers a bit, we can get a simulation of the moon’s orbit!

If you want to play with the code, here it is: Earth Moon Simulation on Codepen.

Collisions & Constraints

While that mathematical model of things I presented does model the physical world, the numerical integration method sadly falls apart in certain situations.

Consider a simulation of a bouncing ball.

The entire state of the simulation is:

Where \( x \) is the ball’s height off the ground and \( v \) is the velocity of the ball. If we drop the ball from a height of 0.8 meters, we have:

If we plot \( x(t) \), it would look something like this:

While the ball is falling, we can construct our derivative function \( f \) fairly easily:

While accelerating only under the influence of gravity, \( a = -g = -9.8 \frac{m}{s^2} \).

But what happens when the ball hits the ground? We know the ball has hit the ground when \( x = 0 \). But in our numerical integration, it’s possible that the ball might be above the ground at one time step, and in the ground the next time step: \( x_i > 0, x_{i+1} < 0 \).

We could solve this by rewinding time to find the time \( t_c \) at which the collision happens \( (t_i < t_c < t_{i+1}) \). But once we’ve found that, we still don’t have a way to define \( \frac{dv}{dt} \) that would result in the velocity instantaneously changing to be upwards instead of downwards.

It’s possible to work this all out by having the collision take a finite amount of time to take place and applying some force \( F \) over that timespan \( \Delta t \), but it’s usually easier to just support some kind of discrete constraints to the simulation.

We can also do more than one iteration of the physics simulation per rendered frame to reduce the amount by which the ball penetrates the floor before bouncing. With that in mind, our core simulation code changes to this:

function runSimulation<T extends Numeric<T>>(
  y0: T,
  f: (t: number, y: T) => T,
  applyConstraints: (y: T) => T,
  iterationsPerFrame: number,
  render: (y: T) => void
) {
    const frameTime = 1 / 60.0
    const h = frameTime / iterationsPerFrame

    function simulationStep(yi: T, ti: number) {
        render(yi)
        requestAnimationFrame(function () {
            for (let i = 0; i < iterationsPerFrame; i++) {
                yi = yi.plus(f(ti, yi).times(h))
                yi = applyConstraints(yi)
                ti = ti + h
            }
            simulationStep(yi, ti)
        })
    }
    simulationStep(y0, 0.0)
}

And then we can implement our bouncing ball like so:

const g = -9.8; // m / s^2
const r = 0.2; // m

class Ball implements Numeric<Ball> {
    constructor(readonly x: number, readonly v: number) { }
    plus(other: Ball) { return new Ball(this.x + other.x, this.v + other.v) }
    times(scalar: number) { return new Ball(this.x * scalar, this.v * scalar) }
}

function f(t: number, y: Ball) {
    const { x, v } = y
    return new Ball(v, g)
}

function applyConstraints(y: Ball): Ball {
    const { x, v } = y
    if (x <= 0 && v < 0) {
        return new Ball(x, -v)
    }
    return y
}

const y0 = new Ball(
    /* x */ 0.8,
    /* v */ 0
)

function render(y: Ball) {
    ctx.clearRect(0, 0, 400, 400)
    ctx.fillStyle = '#EB5757'
    ctx.beginPath()
    ctx.ellipse(200, 400 - ((y.x + r) * 300), r * 300, r * 300, 0, 0, 2 * Math.PI)
    ctx.fill()
}

runSimulation(y0, f, applyConstraints, 30, render)

You can play around with this version of the code here: Bouncing Ball on Codepen.

Implementers beware!

While this general description of simulations has some nice properties, it doesn’t necessarily yield the most high performance simulations. I find it a nice framework to think about the behavior of simulations, though there’s certainly a lot of unnecessary overhead.

The rain simulation that starts this post is, for instance, not implemented using the patterns here. Instead, it was an experiment using the Entity Component System architectural pattern. You can see the source for the rain simulation here: Rain Simulation source on GitHub.

Until next time!

Something about the intersection of math, physics, and programming really strikes a chord with me. Getting simulations up and running and rendering makes for a very special kind of something out of nothing.

SIGGRAPH course notes served as a form of inspiration here, just as they did with the fluid simulation. Check out “An Introduction to Physically Based Modelling” course notes from SIGGRAPH 2001 if you want a much more rigorous look into this material. I have the notes from 1997 linked because it looks like Pixar took the 2001 versions down.

Thanks to Maggie Cai for the beautiful illustration of the couple with the umbrella and for having the patience to meticulously choose colors when I can barely tell the difference between blue and grey.

And in case you were wondering, the illustration was made in Figma. 😃

This post is also available in Japanese: 一から学ぶベジェ曲線.

How do you describe a straight line segment? We might think about a line segment in terms of its endpoints. Let’s call those endpoints \( P_0 \) and \( P_1 \).

To define the line segment rigorously, we might say “the set of all points along the line through \( P_0 \) and \( P_1 \) which lie between \( P_0 \) and \( P_1 \)“, or perhaps this:

$$ L(t) = (1 - t) P_0 + t P_1, 0 \le t \le 1 $$

Conveniently, this definition let’s us easily find the coordinate of the point any portion of the way along that line segment. The midpoint, for instance, lies at \( L(0.5) \).

We can, in fact, linearly interpolate to any value we want between the two points, with arbitrary precision. This allows us to do fancier things, like trace the line by having the \( t \) in \( L(t) \) be a function of time.

If you got this far, you might now be wondering, “What does this have to do with curves?“. Well, it seems quite intuitive that you can precisely describe a line segment with only two points. How might you go about precisely describing this?

It turns out that this particular kind of curve can be described by only 3 points!

This is called a Quadratic Bezier Curve. A line segment, donning a fancier hat, might be called a Linear Bezier Curve. Let’s investigate why.

First, let’s consider what it looks like when we interpolate between \( P_0 \) and \( P_1 \) while simultaneously interpolating between \( P_1 \) and \( P_2 \).

Now let’s linearly interpolate between \( B_{0, 1}(t) \) and \( B_{1, 2}(t) \)…

Notice that the equation for \( B_{0, 1, 2}(t) \) looks remarkably similar to the equations for \( B_{0, 1} \) and \( B_{1, 2} \). Let’s see what happens when we trace the path of \( B_{0, 1, 2}(t) \).

We get our curve!

Higher Order Bezier Curves

Just as we get a quadratic bezier by interpolating between two linear bezier curves, we get a cubic bezier curve by interpolating between two quadratic bezier curves:

You may have a sneaking suspicion at this point that there’s a nice recursive definition lurking here. And indeed there is:

Or, expressed (concisely but inefficiently) in TypeScript, it might look like this:

type Point = [number, number];
function B(P: Point[], t: number): Point {
    if (P.length === 1) return P[0];
    const left: Point = B(P.slice(0, P.length - 1), t);
    const right: Point = B(P.slice(1, P.length), t);
    return [(1 - t) * left[0] + t * right[0],
            (1 - t) * left[1] + t * right[1]];
}
// Evaluate a cubic spline at t=0.7
B([[0.0, 0.0], [0.0, 0.42], [0.58, 1.0], [1.0, 1.0]], 0.7)

Cubic Bezier Curves in Vector Images

As it happens, cubic bezier curves seem to be the right balance between simplicity and accuracy for many purposes. These are the kind of curves you’ll most often see in vector editing tools like Figma.

You can think of the two filled in circles ● as \( P_0 \) and \( P_3 \), and the two diamonds ◇ as \( P_1 \) and \( P_2 \). These are the fundamental building blocks of more complex curved vector constructions.

Font glyphs are specified in terms of bezier curves in TrueType (.ttf) fonts.

The Scalable Vector Graphics (.svg) file format uses bezier curves as one of its two curve primitives, which are used extensively in this:

Cubic Bezier Curves in Animation

While bezier curves have their most obvious uses in representing spacial curves, there’s no reason why they can’t be used to represented curved relationships between other quantities. For instance, rather than relating \( x \) and \(y \), CSS transition timing functions relate a time ratio with an output value ratio.

Cubic bezier curves are one of two ways of expressing timing functions in CSS (steps() being the other). The cubic-bezier(x1, y1, x2, y2) notation for CSS timing functions specifies the coordinates of \( P_1 \) and \( P_2 \) of a cubic bezier curve.

Let’s pretend we’re trying to animate an orange ball moving. In all of these diagrams, the red lines representing time move at a constant speed.

Why Bezier Curves?

Bezier curves are a beautiful abstraction for describing curves. The most commonly used form, cubic bezier curves, reduce the problem of describing and storing a curve down to storing 4 coordinates.

Beyond the efficiency benefits, the effect of moving the 4 control points on the curve shape is intuitive, making them suitable for direct manipulation editors.

Since 2 of the points specify the endpoints of the curve, composing many bezier curves into more complex structures with precision becomes easy. The exact specification of endpoints is always what makes it so convenient in the animation case: the only sensible value of the easing function at \( t = 0\% \) is the initial value, and the only sensible value at \( t = 100\% \) is the final value.

A less obvious benefit is that the line from \( P_0 \) to \( P_1 \) specifies the tangent of the curve leaving \( P_0 \). This means if you have two joint curves with mirrored control points, the slope at the join point is guaranteed to be the same on either side of the join.

A major benefit of mathematical construct like bezier curves is the ability to leverage decades of mathematical research to solve most problems you might run into, completely agnostic to the rest of your problem domain.

For instance, to make this post, I had to learn how to split a bezier curve at a given value of \( t \) in order to animate the curves above. I was quickly able to find a well written article on the matter: A Primer on Bézier Curves: Splitting Curves.

Resources and Further Reading

• A Primer on Bézier Curves in addition to having a description of using deCasteljau’s algorithm to draw and split curves, this free online book seems to be a pretty comprehensive intro.
• Bézier curve on Wikipedia shows many different mathematical formulas of bezier curves beyond the recursive definition shown here. It also contains the original animations that made bezier curves seem so evidently elegant to me.

Also a shoutout to Dudley Storey for his article Make SVG Responsive, which allowed all of the inline SVG in this article to work nicely on mobile.

1. Iteration Speed: A sub-second edit-run cycle
2. Correctness Checking: A compiler that can tell me that my code is probably wrong without me having to run every single codepath
3. Concise syntax: Express intent fluidly without a lot of boilerplate
4. Editing support: Autocompletion, go to definition, go to usage, refactoring support that always works
5. Debugging support: Errors that are easy to debug

First, I want to talk about the tradeoffs that two historical camps make here, then about how these tradeoffs are being bent, and finally how people bend them the opposite way and put themselves in sad, unnecessary hell.

The Dynamically Typed Camp

The dynamically typed camp of languages is where Python, JavaScript, PHP, Ruby, and Scheme live. This is where I’ve spent most of my professional career, including nearly all of the time I spent at Khan Academy.

Iteration Speed

Dynamically typed languages do great here. You edit a file, re-run your program (possibly it automatically restarts itself), and you’re probably back in business in less than a second. This is absolutely wonderful for iterating on ideas quickly and especially for doing things like pixel pushing in UIs.

Correctness Checking

This is where dynamically typed languages fall flat on their faces. Let’s say you have a whole bunch of classes with a method called .render(title), and they all currently take one argument. Now you want to change one of those functions to look like .render(title, description). Someone sends you a code review with this in it, and you ask if they updated all the right places. They respond “I hope so”, and if you have a 100KLOC codebase, with minimal test coverage, then that’s pretty much as good as you’re going to get.

Concise Syntax

Dynamic languages here are generally pretty solid. With no type definitions, there’s just generally less to type. To use opposite extremes here, the difference here is pretty clear.

# Python
a = [1, 2, 3]

// Old versions of Java
ArrayList<Integer> list = new ArrayList<Integer>();
list.add(1);
list.add(2);
list.add(3);

Editing Support

For the same reason as correctness checking, doing refactorings that get more complex than a global string search and replace start to seriously suck in dynamic languages. If you want to delete a method with a common name, or get autocomplete on an object passed into a function, you might be shit out of luck.

You can still get useful autocomplete, but it’ll be at least a little crippled because of the fundamental ambiguity of what a line of code does until it’s actually run in dynamically typed languages.

Debugging Support

The ability to drop down into a REPL with the full power of the language whenever something crashes is actually one of the major reasons I’ve been so happy in the dynamic camp for so long. In Python this looks like pdb.set_trace(), in JavaScript, like debugger, in Ruby, like binding.pry. In the middle of your breakpoint, you can define new functions, invoke arbitrary functions, write data to files – whatever you want.

The Statically Typed Camp

Industry staples like C++, Java, and Objective-C live here with their functional comrades Haskell and OCaml, plus some new company of Go, Swift, Scala, Rust, and Elm. I’m living in this land for the foreseeable at Figma working in C++ and TypeScript.

Iteration Speed

This is arguably the biggest downside of statically typed languages. Type checking, as it turns out, is frequently slow. And since in most of these languages, type resolution is pre-requisite to code generation, slow type checking means slow compiling. Slow compiling means slow iteration time.

While I was an intern at Facebook, I needed to make some changes to WebKit. The compile time on my Macbook was ~15 minutes for every change. Suffice to say, this was not a fun experience. On the upside, I did read most of Pro Git while I was waiting for XCode to build.

Correctness Checking

If your Post and Picture classes both have a .render, and you want to change the signature on Post but not on Picture, you’ve got no troubles in static land. An IDE will make this as easy as right clicking and “Change Signature”. And if you do decide to do it manually because you need to go manually decide on that second argument value at all the new call-sites, no problem – your compiler will quite happily tell you if you done goofed or not.

The level of safety you get here varies wildly by language. Most notable, most compilers don’t work too hard to try to figure out if a pointer is null before telling you your code is A-OK.

In C++, the compiler will quite happily let you do this:

User* a = nullptr;
a->setName("Gertrude")

Haskell and Scala do their best to dodge this problem by not letting you have null, instead representing optional fields explicitly with a Maybe User/Option[User], where it forces you to deal with the fact that it might be missing, and not just assume it’s there.

Editing Support

Statically typed languages kill it here. Since, by definition, the type of every variable must be known without needing to execute the code, your editor can be quite confident which operations are valid on which variables, and helpfully autocomplete them. It can also facilitate things like field renaming, automatic documentation lookup, consistently working go-to definition, and go-to usages.

Debugging Support

My experience varies here, but for the most part have been displeased by my debugging experiences in statically typed languages. While gdb and friends will let you evaluate certain expressions, you lose the ability to do arbitrary manipulations like define debugging helper functions or easily write function invocations on anything templated.

A particularly nasty class of this where you don’t get any interactive console at all to debug is complex compile errors. In C++, clang improved this dramatically, but for code like this:

#include <vector>
#include <algorithm>
int main()
{
    int a;
    std::vector< std::vector <int> > v;
    std::vector< std::vector <int> >::const_iterator it = std::find( v.begin(), v.end(), a );
}

gcc used to output > 15000 characters of errors (see the rest here), which starts like this:

/usr/include/c++/4.6/bits/stl_algo.h: In function ‘_RandomAccessIterator
std::__find(_RandomAccessIterator, _RandomAccessIterator, const _Tp&,
std::random_access_iterator_tag) [with _RandomAccessIterator =
__gnu_cxx::__normal_iterator*, std::vector > >, _Tp = int]’:
/usr/include/c++/4.6/bits/stl_algo.h:4403:45:   instantiated from ‘_IIter std::find(_IIter, _IIter, const _Tp&) [with _IIter = __gnu_cxx::__normal_iterator*, std::vector > >, _Tp = int]’
error_code.cpp:8:89:   instantiated from here
/usr/include/c++/4.6/bits/stl_algo.h:162:4: error: no match for ‘operator==’ in ‘__first.__gnu_cxx::__normal_iterator::operator* [with _Iterator = std::vector*, _Container = std::vector >, __gnu_cxx::__normal_iterator::reference = std::vector&]() == __val’
/usr/include/c++/4.6/bits/stl_algo.h:162:4: note: candidates are:

So that sucks, and I don’t even have anything I can play with to introspect what the issue is.

Stuck

So for an era of programming, it felt like you were kind of stuck between two worlds, each of which had pretty crappy tradeoffs. Then the two camps stopped yelling across the river at each other and started to recognize that the other team was maybe onto something. You see a similar middle ground emerging in the Object Oriented vs. Functional holy war with languages like Scala and Swift taking an OO syntax, functional thinking approach, and JavaScript being kind of accidentally multi-paradigm.

But back to types. Let’s talk about how people are trying to have their cake and eat it too.

The Best of Both Worlds

Type Inference

Something about the following line of Java just feels insulting.

ArrayList<Integer> list = new ArrayList<Integer>();

Why do I need to specify the type information twice? This feels super dumb. More generally, if I do:

a = somefunction()

And the compiler knows the return type of somefunction, I, as the programmer, shouldn’t be forced to tell the computer information it already knows.

The more complete version of this idea is type inference, and the first time I saw it was in Haskell, as concisely explained in Learn You a Haskell for Great Good!: Types and Typeclasses, and is explained with a bit more detail in TypeScript’s Type Inference Documentation.

So we get concise syntax without sacrificing type information.

It’s also now made its way into C++ via the C++11 auto keyword, and is a feature of most modern statically typed languages like Scala, Swift, Rust, and Go.

Decoupling Type Checking from Code Generation

If you define your language very carefully, you can make the compiler output not dependent on the types (i.e. ignore the type information completely), and then run type checking completely separately. The easiest way of defining a language like this is to start with a dynamically typed language and start adding type annotations. Facebook’s Flow does this by adding type annotations to JavaScript.

For instance, a bit of Flow annotated JavaScript might look like this:

// @flow
function bar(x): string {
  return x.length;
}
bar('Hello, world!');

Compilation here is incredibly fast, because all it does is strip the type annotations to produce this:

function bar(x) {
  return x.length;
}
bar('Hello, world!');

While the type checker runs in the background, or possible only on demand. With this type information, you can get better correctness guarantees, and much better IDE support. We get all of this without the normal increase in iteration time that comes with a blocking type checker.

Facebook took a similar approach to type annotating PHP with its language Hack. Python 3.5 introduced progressive typing too, as described in the typing module.

An interesting side effect of having a compile target that closely resembles the source is that you get all the benefits of the interactive debugging mentioned before. Microsoft’s TypeScript takes a very similar approach to Flow, except that you can set a flag to localize type checking to file-internal type checking, making the assumption that all imports are of the correct type, which speeds up type checking considerably.

Linting

Linting is a form of static analysis that happens outside of a compiler, and typically on a dynamic language. One of the earliest ones I’d heard of was Douglas Crockford’s JSLint that, despite the possibly dynamic nature of your program, might be able to confidently point out mistakes. There are countless tools that do this for various languages, and I go into more depth about the value of them in Linters as Invariants. This gives you a small subset of the correctness guarantees that you get from a statically typed language, like the guarantee that you aren’t using a variable that isn’t declared anywhere, but typically isn’t very helpful for inter-file analysis.

Good linters, like good IDE support, will allow for automatic fixing of errors, like with the ESLint --fix flag.

Faster Compilers

This one is pretty self explanatory. If your compiler is super fast, the edit-run iteration time isn’t an issue. Boom.

To do this properly, you need to design the language carefully with that as a design goal. Go did this.

Go is an attempt to combine the ease of programming of an interpreted, dynamically typed language with the efficiency and safety of a statically typed, compiled language. It also aims to be modern, with support for networked and multicore computing. Finally, working with Go is intended to be fast: it should take at most a few seconds to build a large executable on a single computer.

– Go FAQ

Better Compiler Error Messages

There have been numerous attempts to make debugging compilation errors a non-issue by having sensible human-readable error messages, notably in Elm, which Evan Czaplicki describes in his post Compiler Errors for Humans. The cool thing about this is that it’s seeing adoption by slightly more main-stream languages like Rust, as Jonathan Turner explains in Shape of errors to come. Much earlier, improvements to C++ error messages were a selling point of clang over gcc as described in Expressive Diagnostics. But Elm and Rust take it steps further.

So if you have a smart enough compiler, ideally one that can suggest to you how to fix your problem, then some of the frustration of debugging those issues melt away.

All the Downside

Conversely, some very old techniques in the static world, and some very new techniques in the dynamic world put you in a world of pain by taking downsides from both camps.

Runtime Failures in a Static Language

One of my least pleasant recurring memories working as an intern at Square in Java was waiting a few minutes for a Java build, only to have it crash on boot with a runtime error. This was happening because of fancy dynamic runtime dependency injection with Guice, which is why Square later wrote a library that does it at compile time called Dagger to fix the problem.

The more general Bad Situation to avoid here are things that pass type checking, cause runtime crashes, and require a recompile to fix.

Examples of this include things like:

• using runtime dependency injection like Guice
• using raw pointers instead of smart pointers in C++.
• liberal use of void* in C/C++ or using the Object type in Java then doing run-time downcasts all over the place

Now you’re waiting for 15s for compilation, only to find that your code crashes on boot repeatedly.

Slow Transpiling/Bundling in a Dynamic Language

Just as the example above takes a static language that should be safe and makes it unsafe, you can take a dynamic language and make it slow to iterate on! The most common culprit of this is doing complex transpilation for a lot of code, and doing that on every code change.

“Transpiling” started gaining momentum in the web development world when compile-to-css languages like Less and Sass, and compile-to-js languages like CoffeeScript rolled around, but really blew up when React, webpack and babel started becoming a trio of choice.

The idea of using a more expressive languages than CSS and JavaScript to write safer, more readable code is wonderful, but if you’re not careful, you’ve now managed to inherit the increased edit-run cycle time of a statically typed language without inheriting any of the correctness guarantees.

Congratulations, you now have an unresponsive, unmaintainable mess.

Closing thoughts

Overall, I’m pretty happy with the direction that things are going in PL world. I hope much of the near future will be built on TypeScript (with --strictNullChecks) instead of JavaScript on the front end, Rust instead of C++, Scala instead of Java, Go instead of Node.js, and Swift instead of Objective-C.

I haven’t had the chance to play much with Go, Rust, or Swift, but things sound kinda rosy.

If you like thinking about language tradeoffs and want stories more informed by experience, you should read through Steve Yegge’s Is Weak Typing Strong Enough?.

EDIT: My fellow uWaterloo Software Engineering 2014 classmate, Ming-Ho Yee, is a PhD student in programming language design, so he naturally had some things to say about this post. You can read his thoughts in his response.

Note: The demos in this post rely on WebGL features that might not be implemented in mobile browsers.

About a year and a half ago, I had a passing interest in trying to figure out how to make a fluid simulation. At the time, it felt just a bit out of my reach, requiring knowledge of shaders, vector calculus, and numerical computation that were all just a little bit past my grasp. At the time, I was working through the Fluid Simulation Course Notes from SIGGRAPH 2007, and was struggling with the math. Now armed with a bit more knowledge and a lot more time, and with the help of other less dense resources like GPU Gems Chapter 38. Fast Fluid Dynamics Simulation on the GPU, I was finally able to figure out enough to get something working. I am still a beginner at simulations like this, and I’m going to brazenly ignore things like numerical stability, but hopefully I can help leapfrog you past a few places I got stuck.

We’re going to work with the simplest 2D fluid simulation, where the entire area is full of fluid, and we’re going to ignore viscosity.

The Velocity Field

As compared to a rigid, unrotating solid, where every bit of the thing has to be moving in the same direction at the same speed, each bit of a fluid might be moving differently. One way to model this is to use a vector field representing velocity. For any given \( (x, y) \) coordinate, this field will tell you the velocity of the fluid at that point.

$$ \vec u(x, y) = (u_x, u_y) $$

A nice way to get an intuition about what a given field looks like is to sample the function in a grid of points, then draw arrows starting at each grid point whose size and orientation are dictated by the value of the function at that point. For the purposes of this post, we’re always going to be working over the domain \( x \in [-1, 1] \), \( y \in [-1, 1] \).

For instance, here’s a very simple field \( \vec u(x, y) = (1, 0) \) representing everything moving at a constant speed to the right.

And here’s a more interesting one \( \vec u(x, y) = (x, y) \) where things move away from the origin, increasing in speed the farther away from the origin they are.

We’re going to play with this one, \( \vec u(x, y) = \left( \sin (2 \pi y), \sin (2 \pi x) \right) \), since it creates some interesting visual results once we start making the fluid move accordingly.

For a more thorough introduction to vector fields, check out the Introduction to Vector Fields video on Khan Academy. The rest of the videos on multivariate calculus might prove useful for understanding concepts in fluid flow too.

Now then, let’s get things moving.

Advection

Advection is the transfer of a property from one place to another due to the motion of the fluid. If you’ve got some black dye in some water, and the water is moving to the right, then surprise surprise, the black dye moves right.

If the fluid is moving in a more complex manner, that black dye will get pulled through the liquid in a more complex manner.

Before we dive into how advection works, we need to talk a bit about the format of the data underlying these simulations.

The simulation consist of two fields: color and velocity. Each field is represented by a two dimensional grid. For simplicity, we use the same dimensions as the output pixel grid.

Previously, I described the velocity field as an analytical function \( \vec u(x, y) \). In practice, that analytical function is only used to initialize the grid values.

The simulation runs by stepping forward bit-by-bit in time, with the state of the color and velocity grids depending only on the state of the color and velocity grids from the previous time step. We’ll use \( \vec u(\vec p, t) \) to represent the velocity grid at 2d position \( \vec p \) and time \( t \), and \( \vec c(\vec p, t) \) to represent the color in the same manner.

So how do we move forward in time? Let’s just talk about how the color field changes for now. If we consider each grid point as a little particle in the fluid, then one approach is to update the color of the fluid where that particle will be, one time step in the future.

In order to run these simulations in real-time at high resolution, we want to implement them on the GPU. It turns out that this method of updating the value at the new location of the particle is difficult to implement on the GPU.

First, the position we want to write, \( \vec p + \vec u(\vec p, t) \Delta t \) might not lie on a grid point, so we’d have to distribute the impact of the write across the surrounding grid points. Second, many of our imaginary particles might end up in the same place, meaning we need to analyze the entire grid before we decide what the new values of each grid point might be.

So, instead of figuring out where our imaginary particles at the grid points go to, we’ll figure out where they came from in order to calculate the next time step.

With this scheme, we only need to write to a single grid point, and we don’t need to consider the contributions of imaginary particles coming from multiple different places.

The last teensy hurdle is figuring out the value of \( \vec c(\vec p - \vec u(\vec p, t) \Delta t, t ) \), since \( \vec p - \vec u(\vec p, t) \) might not be at a grid point. We can hop this hurdle using bilinear interpolation on the surrounding 4 grid points (the ones linked by the dashed grey rectangle above).

Advecting the Velocity Field

Barring a bizarre sequence of perfectly aligned fans underneath the liquid, there’s no reason why the velocity field wouldn’t change over time. Just as black ink would move through the fluid, so too will the velocity field itself! Just as we can advect \( \vec c \) through \( \vec u \), we can also advect \( \vec u \) through itself!

Intuitively you can think of it this way: a particle moving in a certain direction will continue moving in that direction, even after it’s moved.

Since we’re storing velocity in a grid just like we did with color, we can use the exact same routine to advect velocity through itself. Below, watch the velocity change over time, with an initial velocity field of \( \vec u = (1, \sin(2 \pi y)) \).

If you tried playing around with this, and saw a bunch of weird hard edges and might’ve thought to yourself “I don’t think fluids work like that…”, you’d be right. We’re missing an important ingredient, but before we look at the solution, let’s take a closer look at the problem.

Divergent Fields

Something about the velocity field below makes this intuitively not feel like a fluid. Fluids just don’t behave like this.

Same problem with this one…

If you look at where the arrows are pointing in each of the above 2 simulations, you’ll see that there are spots where the all the arrows point away from that spot, and others where all the arrows point toward that spot. Assuming the volume of the liquid is staying constant, the density of the fluid has to be changing for such a velocity field to be possible.

Water is roughly incompressible. That means that at every spot, you have to have the same amount of fluid entering that spot as leaving it.

Mathematically, we can represent this fact by saying a field is divergence-free. The divergence of a velocity field \( \vec u \), indicated with \( div(\vec u) \) or \( \nabla \cdot \vec u \), is a measure of how much net stuff is entering or leaving a given spot in the field. For our 2D velocity field, it’s defined like this:

The first of the two not-very-fluidy fields above has an equation \( \vec u(x, y) = (x, y) \). Taking the divergence, we find:

This positive value tells us that, in all places, more stuff is leaving that point than entering it. In physical terms, this means that the density is decreasing uniformly everywhere.

The other not-very-fluidy field has an equation \( \vec u(x, y) = \sin(2 \pi x), 0) \). If we look at its divergence, we see:

Which tells us that in some places, density is increasing (where \( \nabla \cdot \vec u < 0 \)), and in others, density is decreasing (where \( \nabla \cdot \vec u > 0 \)).

Doing the same operation on the more fluidy looking swirly velocity field \( \vec u = (\sin ( 2 \pi y), \sin ( 2 \pi x ) \) that you saw in the section about advection, we discover \( \nabla \cdot \vec u = 0 \).

An incompressible fluid will have a divergence of zero everywhere. So, if we want our simulated fluid to look kind of like a real fluid, we better make sure it’s divergence-free.

Since our velocity field undergoes advection and can be influenced by clicking and dragging around the fluid, having an initially divergence-free velocity field isn’t enough to guarantee that the field will continue to be divergence-free. For example, if we take our swirly simulation and start advecting the velocity field through itself, we end up with something divergent:

So we need a way of taking a divergent field and making it divergence-free. To understand what force makes that happen in the real world, we need to talk about some honest-to-goodness physics.

Navier-Stokes

The Navier-Stokes equations describe the motion of fluids. Here are the Navier-Stokes equations for incompressible fluid flow:

Where \( \vec u \) is the velocity field, \( \rho \) is density, \( p \) is pressure, \( \nu \) is the kinematic viscosity, and \( \vec F \) is external forces acting upon the fluid.

Since we’re pretending the viscosity of our fluid is zero, we can drop the \( \nu \) term in the first equation. In our simple simulation, external forces are only applied by dragging the mouse, so we’ll ignore that term for now, opting to allow it to influence the velocity field directly.

Dropping those terms, we’re left with the following:

We can expand this to its partial derivative form, expanding vector components to leave us with only scalar variables.

Remembering that these fields are all functions on \( (x, y, t) \), we can approximate the partial derivatives with finite differences. For instance, we can approximate the partial derivative of \( u_x \) with respect to \( t \) like so:

Because the procedure ends up being the same for both components, we’ll focus on only the \( x \) component here. Applying finite differences to all of the partial derivatives, we have this:

Ultimately what we want is \( \vec u(x, y, t + \Delta t) \), which will tell us, for a given point, what the velocity will be at the next time step. So let’s solve for that by rearranging the big long formula above:

If you look at the first three terms in this expression, what does it look like they conceptually represent? It looks like they represent the next velocity after we’ve taken into account changes due to the motion of the fluid itself. That sounds an awful like advection as discussed earlier. In fact, it will work quite well if we substitute the velocity field after it’s undergone advection. We’ll call the advected velocity field \( \vec u ^ a \). So now we have:

So after all of that, we have an equation that relates the velocity field at the next time tick to the current velocity field after it’s undergone advection, followed by application of pressure.

We know that a divergence-free field that undergoes advection isn’t necessarily still divergence-free, and yet we know that the Navier-Stokes equations for impressible flow represent divergence-free velocity fields, so therefore we have our answer about what in nature prevents the velocity field from becoming divergent: pressure!

Solving for Pressure

Now that we have an equation that relates \( \vec u \) to \( p \). This is where the math gets messy. We start from the second Navier-Stokes equation for incompressible flow, applied at time \( t + \Delta t \), and apply finite differences again:

Here, we can substitute our equations for \( \vec u \) expressed in terms of \( \vec u ^ a \) and \( p \) to get this monster:

Rearranging to have all of the \( p \) terms on the left and all the \( \vec u ^ a \) terms on the right, and multiplying both sides by \( 2 \epsilon \), we have:

Note: the above expression is a scalar expression, despite being laid out in a somewhat vector-y form.

At this point, it’s helpful to remember that, for the purposes of the simulation, we’re not interested in knowing the value of \( p \) everywhere: we only care about knowing its value at enough places to calculate the value of the velocities at the grid points. To meet that end, we can similarly calculate \( p \) on the grid. To accomplish this, we can make \( \epsilon \) the distance between adjacent grid cells.

The above equation yields a new equation for every \( (x, y) \) of a grid point we substitute. For the purposes of discussion, let’s assume that the gap between adjacent cells is 0.1 units, so \( \epsilon = 0.1 \). Let’s examine what the equation yields for \( (x, y) = (0.3, 0.7) \).

All the values on the right hand side of this equation are known, and on the left we have 5 unknowns: the value of \( p \) at 5 different grid locations.

If we repeat this process and evaluate \( (x, y) \) at every grid point, we get one equation with 5 unknowns for each grid location. If our grid has \( n \times m \) grid locations in it, then we have \( n \times m \) equations, each with 5 unknowns.

If you’re wondering about what’s happening at the edges, we’re going to lazily side-step that question by making our grid wrap around: if you ask for the velocity past the bottom edge, you’ll get a value near the top edge.

Before we move on, our notation is getting a bit clunky, so let’s clean it up a tad since we know we’re working on a grid. For the next part, we’ll say \( p_{i,j} = p(i \epsilon, j \epsilon, t) \), and we’ll stick all the known values together into a value \( d \) (for divergence), like so:

With this nicer notation, we can express the system of equations on pressure that we’re trying to solve like so:

Iteratively Solving the Pressure Equation

Solving for \( p_{i, j} \) for every grid point analytically would be an enormous mess. Instead, we’re going to use an iterative method of solving this system of equations, where each iteration provides values closer and closer to a real solution. We’re going to use the Jacobi Method.

In the Jacobi method, we first rearrange our equation to isolate one term, like so:

Next, we make an initial guess for all of our unknowns. We’ll call this initial guess \( p_{i,j}^{(0)} \), and just set it to 0 everywhere.

Here’s where the iteration comes in: our next guess, \( p_{i, j}^{(1)} \) is obtained by plugging in our initial guess into the above formula:

And, more generally, each iteration relies upon the previous one:

You would usually run this until the values of one iteration are equal to the values from the previous iteration, rounded to a certain accuracy. For our purposes, we’re more interested in this running in a consistent period of time, so we’ll arbitrarily run this for 10 iterations, and hope the result is accurate enough to look realistic.

For a bit of intuition on why this converges to a solution, check out Algebraic Pavel’s answer on Math Exchange.

All Together Now!

Phew! That was a lot to get through. Now let’s put it all together. Roughly, as pseudo-code, here’s our whole simulation:

initialize color field, c
initialize velocity field, u

while(true):
    u_a := advect field u through itself
    d := calculate divergence of u_a
    p := calculate pressure based on d, using jacobi iteration
    u := u_a - gradient of p
    c := advect field c through velocity field u
    draw c
    wait a bit

Here are the key formulas for those steps on grid coordinates \( (i, j) \), uncluttered by derivations:

Advecting field \( vec u \) through itself:

Divergence of \( \vec u_a \) (multiplied by constant terms):

Pressure calculation Jacobi iteration step, with \( p_{i, j}^{(0)} = 0 \):

Subtracting the pressure gradient from the advected velocity field:

Advecting the color field through the final velocity field:

Implementation

Pulling all those steps together, you can make something like this! Woohoo! When I got this working for the first time, I was pretty ecstatic.

I won’t delve too far into the implementation, but you can have a look at it yourself: fluid-sim.js. It relies upon the elegant lightgl.js, which is an abstraction layer on top of WebGL that makes it much nicer to work with. Unlike THREE.js, it doesn’t make any assumptions about you wanting any concept of a camera or lighting or that you’re working in 3D at all.

The key technique for running the simulation efficiently is doing all the hard work on the GPU. To meet this need, all of the computations are done via the render to texture technique, ping-ponging which texture is being rendered to facilitate reading and writing to the same conceptual texture (e.g. reading from the velocity field and writing to the velocity field representing the next time step).

Each one of the major components of the algorithm is implemented in a separate shader. There’s a shader for advection, a shader for calculating the divergence, one for a single iteration of the Jacobi method, and another for subtracting the pressure gradient from the advected velocity.

References

To make this, I had to draw from a lot of difference references, many of which are linked inline in the post.

• Fluid Simulation Course Notes from SIGGRAPH 2007: Now a textbook, this is a pretty mathematically dense tutorial. It took me 4 or 5 times reading through most sections to make sense of it, and ultimately I only understood parts of it after I did the derivations myself. It uses the more complex conjugate gradient method instead of the Jacobi method to solve the system of pressure equations, which I got completely lost in, and abandoned. It delves into a lot of arguments about numerical accuracy and uses a more complex grid layout than I did, which I still don’t follow fully. It also has resources for other kinds of fluid simulations, like heightfield simulation, and smoothed particle hydrodynamics.

• GPU Gems Chapter 38. Fast Fluid Dynamics Simulation on the GPU: This was the single most useful reference I found, and describes something very similar to this post. It walks through specific implementation ideas, and gave me a much better intuition for advection. Some of the math (or at least the notation) seems shaky here. I think the Gaussian “splat” formula is missing a negative sign inside of the \( exp() \), and I’m not sure what the notation \( (\vec u \cdot \nabla) u_x \) means in the first Navier-Stokes equation, since \( \nabla \cdot \vec u = 0 \) in the second equation.

• “Elementary Linear Algebra” by Ron Larson, Section 10.2: Iterative Methods for Solving Linear Systems. This had a much clearer explanation of the Jacobi method than the GPU Gems chapter that allowed me to derive the pressure solve iteration myself. The full textbook can be found here: “Elementary Linear Algebra” on amazon.com.

• Jonas Wagner’s fluid simulation on canvas, and particularly the source for it (fluid.js) were helpful for understanding what a full solution actually looks like. It’s also how I found the GPU Gems article in the first place. Jonas went on later to reimplement his solution in WebGL: WebGL Fluid Simulation.

I’ve always been fascinated by demoscenes: short, real-time generated audiovisual demos, usually in very, very small executables. This one, by an artist named “reptile”, comes from a 4KB compiled executable. No external assets (images, sound clips, etc.) are used – it’s all in that 4KB.

To get a intuitive grip on how small 4KB is, a 1080p video that this produces is 40MB, or ~10,000 times larger than the executable that produces it. Keep in mind that the executable also contains the code to generate the music.

One of the techniques used in many demo scenes is called ray marching. This algorithm, used in combination with a special kind of function called “signed distance functions”, can create some pretty damn cool things in real time.

Signed Distance Functions

Signed distance functions, or SDFs for short, when passed the coordinates of a point in space, return the shortest distance between that point and some surface. The sign of the return value indicates whether the point is inside that surface or outside (hence signed distance function). Let’s look at an example.

Consider a sphere centered at the origin. Points inside the sphere will have a distance from the origin less than the radius, points on the sphere will have distance equal to the radius, and points outside the sphere will have distances greater than the radius.

So our first SDF, for a sphere centered at the origin with radius 1, looks like this:

Let’s try some points:

Great, \( (1, 0, 0) \) is on the surface, \( (0, 0, 0.5) \) is inside the surface, with the closest point on the surface 0.5 units away, and \( (0, 3, 0) \) is outside the surface with the closest point on the surface 2 units away.

When we’re working in GLSL shader code, formulas like this will be vectorized.
Using the Euclidean norm, the above SDF looks like this:

Which, in GLSL, translates to this:

float sphereSDF(vec3 p) {
    return length(p) - 1.0;
}

For a bunch of other handy SDFs, check out Modeling with Distance Functions.

The Raymarching Algorithm

Once we have something modeled as an SDF, how do we render it? This is where the ray marching algorithm comes in!

Just as in raytracing, we select a position for the camera, put a grid in front of it, send rays from the camera through each point in the grid, with each grid point corresponding to a pixel in the output image.

The difference comes in how the scene is defined, which in turn changes our options for finding the intersection between the view ray and the scene.

In raytracing, the scene is typically defined in terms of explicit geometry: triangles, spheres, etc. To find the intersection between the view ray and the scene, we do a series of geometric intersection tests: where does this ray intersect with this triangle, if at all? What about this one? What about this sphere?

Aside: For a tutorial on ray tracing, check out scratchapixel.com. If you’ve never seen ray tracing before, the rest of this article might be a bit tricky.

In raymarching, the entire scene is defined in terms of a signed distance function. To find the intersection between the view ray and the scene, we start at the camera, and move a point along the view ray, bit by bit. At each step, we ask “Is this point inside the scene surface?”, or alternately phrased, “Does the SDF evaluate to a negative number at this point?“. If it does, we’re done! We hit something. If it’s not, we keep going up to some maximum number of steps along the ray.

We could just step along a very small increment of the view ray every time, but we can do much better than this (both in terms of speed and in terms of accuracy) using “sphere tracing”. Instead of taking a tiny step, we take the maximum step we know is safe without going through the surface: we step by the distance to the surface, which the SDF provides us!

In this diagram, \( p_0 \) is the camera. The blue line lies along the ray direction cast from the camera through the view plane. The first step taken is quite large: it steps by the shortest distance to the surface. Since the point on the surface closest to \( p_0 \) doesn’t lie along the view ray, we keep stepping until we eventually get to the surface, at \( p_4 \).

Implemented in GLSL, this ray marching algorithm looks like this:

float depth = start;
for (int i = 0; i < MAX_MARCHING_STEPS; i++) {
    float dist = sceneSDF(eye + depth * viewRayDirection);
    if (dist < EPSILON) {
        // We're inside the scene surface!
        return depth;
    }
    // Move along the view ray
    depth += dist;

    if (depth >= end) {
        // Gone too far; give up
        return end;
    }
}
return end;

Combining that with a bit of code to select the view ray direction appropriately, the sphere SDF, and making any part of the surface that gets hit red, we end up with this:

Voila, we have a sphere! (Trust me, it’s a sphere, it just has no shading yet.)

For all of the example code, I’m using http://shadertoy.com/. Shadertoy is a tool that lets you prototype shaders without needing to write any OpenGL/WebGL boilerplate. You don’t even write a vertex shader – you just write a fragment shader and watch it go.

The code is commented, so you should go check it out and experiment with it. At the top of the shader code for each part, I’ve left some challenges for you to try to test your understanding. To get to the code, hover over the image above, and click on the title.

Surface Normals and Lighting

Most lighting models in computer graphics use some concept of surface normals to calculate what color a material should be at a given point on the surface. When surfaces are defined by explicit geometry, like polygons, the normals are usually specified for each vertex, and the normal at any given point on a face can be found by interpolating the surrounding vertex normals.

So how do we find surface normals for a scene defined by a signed distance function? We take the gradient! Conceptually, the gradient of a function \( f \) at point \( (x, y, z) \) tells you what direction to move in from \( (x, y, z) \) to most rapidly increase the value of \( f \). This will be our surface normal.

Here’s the intuition: for a point on the surface, \( f \) (our SDF), evaluates to zero. On the inside of that surface, \( f \) goes negative, and on the outside, it goes positive. So the direction at the surface which will bring you from negative to positive most rapidly will be orthogonal to the surface.

The gradient of \( f(x, y, z) \) is written as \( \nabla f \). You can calculate it like so:

But no need to break out the calculus chops here. Instead of taking the real derivative of the function, we’ll do an approximation by sampling points around the point on the surface, much like how you learned to calculate slope in a function as rise-over-run before you learned how to do derivatives.

/**
 * Using the gradient of the SDF, estimate the normal on the surface at point p.
 */
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        sceneSDF(vec3(p.x + EPSILON, p.y, p.z)) - sceneSDF(vec3(p.x - EPSILON, p.y, p.z)),
        sceneSDF(vec3(p.x, p.y + EPSILON, p.z)) - sceneSDF(vec3(p.x, p.y - EPSILON, p.z)),
        sceneSDF(vec3(p.x, p.y, p.z  + EPSILON)) - sceneSDF(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

Armed with this knowledge, we can calculate the normal at any point on the surface, and use that to apply lighting with the Phong reflection model from two lights, and we get this:

By default, all of the animated shaders in this post are paused to prevent it from making your computer sound like a jet taking off. Hover over the shader and hit play to see any animated effects.

Moving the Camera

I won’t dwell on this too long, because this solution isn’t unique to ray marching. Just as in raytracing, for transformations on the camera, you transform the view ray via transformation matrices to position and rotate the camera. If that doesn’t mean anything to you, you want to follow through the ray tracing tutorial on scratchapixel.com, or perhaps check out this blog post on codinglabs.net.

Figuring out how to orient the camera based on a series of translations and rotations isn’t always terribly intuitive though. A much nicer way to think about it is “I want the camera to be at this point, looking at this other point.” This is exactly what gluLookAt is for in OpenGL.

Inside a shader, we can’t use that function, but we can look at the man page by running man gluLookAt, take a peek at how it calculates its own transformation matrix, and then make our own in GLSL.

/**
 * Return a transformation matrix that will transform a ray from view space
 * to world coordinates, given the eye point, the camera target, and an up vector.
 *
 * This assumes that the center of the camera is aligned with the negative z axis in
 * view space when calculating the ray marching direction.
 */
mat4 viewMatrix(vec3 eye, vec3 center, vec3 up) {
	vec3 f = normalize(center - eye);
	vec3 s = normalize(cross(f, up));
	vec3 u = cross(s, f);
	return mat4(
		vec4(s, 0.0),
		vec4(u, 0.0),
		vec4(-f, 0.0),
		vec4(0.0, 0.0, 0.0, 1)
	);
}

Since spheres look the same from all angles, I’m switching to a cube here. Placing the camera at \( (8, 5, 7) \) and pointing it at the origin using our new viewMatrix function, we now have this:

Constructive Solid Geometry

Constructive solid geometry, or CSG for short, is a method of creating complex geometric shapes from simple ones via boolean operations. This diagram from WikiPedia shows what’s possible with the technique:

CSG is built on 3 primitive operations: intersection ( \( \cap \) ), union ( \( \cup \) ), and difference ( \( - \) ).

It turns out these operations are all concisely expressible when combining two surfaces expressed as SDFs.

float intersectSDF(float distA, float distB) {
    return max(distA, distB);
}

float unionSDF(float distA, float distB) {
    return min(distA, distB);
}

float differenceSDF(float distA, float distB) {
    return max(distA, -distB);
}

If you set up a scene like this:

float sceneSDF(vec3 samplePoint) {
    float sphereDist = sphereSDF(samplePoint / 1.2) * 1.2;
    float cubeDist = cubeSDF(samplePoint) * 1.2;
    return intersectSDF(cubeDist, sphereDist);
}

Then you get something like this (see section below about scaling to see where the division and multiplication by 1.2 comes from).

In this same Shadertoy, you can play around with the union and difference operations too if you edit the code.

It’s interesting to consider the SDF produced by these binary operations to try to build an intuition for why they work.

Remember that the region where an SDF is negative represents the interior of the surface. For the above intersection, the \( sceneSDF \) can only be negative if both \( cube(\vec p) \) and \( sphere(\vec p) \) are negative. Which means that we only consider a point inside the scene surface if it’s inside both the cube and the sphere. Which is exactly the definition of CSG intersection!

The same kind of logic applies to union. If either function is negative, the resulting scene SDF will be negative, and therefore inside the surface.

The difference operation was the trickiest for me to wrap my head around.

What does the negation of an SDF mean?

If you think again about what the negative and positive region of the SDF mean, you can see that the negative of an SDF is an inversion of the inside and outside of a surface. Everything that was considered inside the surface is now considered outside and vice versa.

This means you can consider the difference to be the intersection of the first SDF and the inversion of the second SDF. So the resulting scene SDF is only negative when the first SDF is negative and the second SDF is positive.

Switching back to geometric terms, that means that we’re inside the scene surface if and only if we’re inside the first surface and outside the second surface – exactly the definition of CSG difference!

Model Transformations

Being able to move the camera gives us some flexibility, but being able to move individual parts of the scene independently of one another certainly gives a lot more. Let’s explore how to do that.

Rotation and Translation

To translate or rotate a surface modeled as an SDF, you can apply the inverse transformation to the point before evaluating the SDF.

Just as you can apply different transformations to different meshes, you can apply different transformations to different parts of the SDF – just only send the transformed ray to the part of the SDF you’re interested in. For instance, to make the cube bob up and down, leaving the sphere in place, but still taking the intersection, you can do this:

float sceneSDF(vec3 samplePoint) {
    float sphereDist = sphereSDF(samplePoint / 1.2) * 1.2;
    float cubeDist = cubeSDF(samplePoint + vec3(0.0, sin(iGlobalTime), 0.0));
    return intersectSDF(cubeDist, sphereDist);
}

Shadertoy reference note: iGlobalTime is a uniform variable set by Shadertoy that contains the number of seconds since playback started.

If you do transformations like this, is the resulting function still a signed distance field? For rotation and translation, it is, because they are “rigid body transformations”, meaning they preserve the distances between points.

More generally, you can apply any rigid body transformation by multiplying the sampled point by the inverse of your transformation matrix.

For instance, if I wanted to apply a rotation matrix, I could do this:

mat4 rotateY(float theta) {
    float c = cos(theta);
    float s = sin(theta);

    return mat4(
        vec4(c, 0, s, 0),
        vec4(0, 1, 0, 0),
        vec4(-s, 0, c, 0),
        vec4(0, 0, 0, 1)
    );
}

float sceneSDF(vec3 samplePoint) {
    float sphereDist = sphereSDF(samplePoint / 1.2) * 1.2;

    vec3 cubePoint = (invert(rotateY(iGlobalTime)) * vec4(samplePoint, 1.0)).xyz;

    float cubeDist = cubeSDF(cubePoint);
    return intersectSDF(cubeDist, sphereDist);
}

…but if you’re using WebGL here, there’s no built in matrix inversion routine in GLSL for now, but you can just do the opposite transform. So the above scene function changes to the equivalent:

float sceneSDF(vec3 samplePoint) {
    float sphereDist = sphereSDF(samplePoint / 1.2) * 1.2;

    vec3 cubePoint = (rotateY(-iGlobalTime) * vec4(samplePoint, 1.0)).xyz;

    float cubeDist = cubeSDF(cubePoint);
    return intersectSDF(cubeDist, sphereDist);
}

For more transformation matrices, refer to any intro to graphics textbook, or check out these slides: 3D Affine transforms.

Uniform Scaling

Okay, let’s get back to this weird scaling trick we glossed over before:

float sphereDist = sphereSDF(samplePoint / 1.2) * 1.2;

The division by 1.2 is scaling the sphere up by 1.2x (remember that we apply the inverse transform to the point before sending it to the SDF). But why do we multiply by the scaling factor afterwards? Let’s examine doubling the size for the sake of simplicity.

float sphereDist = sphereSDF(samplePoint / 2) * 2;

Scaling is not a rigid body transformation – it doesn’t preserve the distance between points. If we transform \( (0, 0, 1) \) and \( (0, 0, 2) \) by dividing them by 2 (which results in a uniform upscaling of the model), then the distance between the points switches from 1 to 0.5.

So when we sample our scaled point in sphereSDF, we’d end up getting back half the distance that the point really is from the surface of our transformed sphere. The multiplication at the end is to compensate for this distortion.

Interestingly, if we try this out in a shader, and use no scale correction, or a smaller value for scale correction, the exact same thing is rendered. Why?

// All of the following result in an equivalent image
float sphereDist = sphereSDF(samplePoint / 2) * 2;
float sphereDist = sphereSDF(samplePoint / 2);
float sphereDist = sphereSDF(samplePoint / 2) * 0.5;

Notice that regardless of how we scale the SDF, the sign of the distance returned stays the same. The sign part of “signed distance field” is still working, but the distance part is now lying.

To see why this is a problem, we need to re-examine how the ray marching algorithm works.

Recall that at every step of the ray marching algorithm, we want to move a distance along the view ray equal to the shortest distance to the surface. We predict that shortest distance using the SDF. For the algorithm to be fast, we want the those steps to be as large as possible, but if we undershoot, the algorithm still works, it just requires more iterations.

But if we overestimate distance, we have a real problem. If we try to scale down the model without correction, like this:

float sphereDist = sphereSDF(samplePoint / 0.5);

Then the sphere disappears completely. If we overestimate the distance, our raymarching algorithm might step past the surface, never finding it.

For any SDF, we can safely scale it uniformly like so:

float dist = someSDF(samplePoint / scalingFactor) * scalingFactor;

Non-uniform scaling and beyond

If we want to scale a model non-uniformly, how can we safely avoid the distance overestimation problem described in the scaling section above? Unlike in uniform scaling, we can’t exactly compensate for the distance distortion caused by the transform. It was possible in uniform scaling because all dimensions were scaled equally, so regardless of where the closest point on the surface to the sampling point is, the scaling compensation will be the same.

But for non-uniform scaling, we need to know where the closest point on the surface is to know how much to correct the distance by.

To see why this is the case, consider the SDF for the unit sphere, scaled to half its size along the x axis, with its other dimensions preserved.

If we evaluate the SDF at \( (0, 2, 0) \), we get back a distance of 1 unit. This is correct: the closest point on the surface of the sphere is \( (0, 1, 0) \). But if evaluate at \( (2, 0, 0) \), we get back a distance of 3 units, which is not right. The closest point on the surface is \( (0.5, 0, 0) \), yielding a world-coordinate distance of 1.5 units.

So, just as in uniform scaling, we need to correct the distance returned by the SDF to avoid overestimating the distance, but by how much? The overestimation factor varies depending on where the point is and where the surface is.

Since it’s usually okay to underestimate the distance, we can just multiply by the smallest scaling factor, like so:

float dist = someSDF(samplePoint / vec3(s_x, s_y, s_z)) * min(s_x, min(s_y, s_z));

The principle for other non-rigid transformations is the same: as long as the sign is preserved by the transformation, you just need to figure out some compensation factor to ensure that you’re never overestimating the distance to the surface.

Putting it all together

With the primitives in this post, you can now create some pretty interesting, complex scenes. Combining those with a simple trick of using the normal vector as the ambient/diffuse component of the material, and you can create something like the shader at the start of the post. Here it is again.

References

There is way more to learn about rendering signed distance functions. One of the most profilic writers on the subject is Inigo Quilez. I learned most of the content of this post from reading his website and his shader code. He’s also one of the co-creators of Shadertoy.

Some of the interesting SDF related material from his site that I didn’t cover at all includes smooth blending between surfaces and soft shadows.

Other references:

• GPU Gems 2: Chapter 8. Per-Pixel Displacement Mapping with Distance Functions
• Raymarching Sample Code on Shadertoy
• Modeling with Distance Functions