Jakub Konka

Disable ASLR when debugging with LLDB on macOS

2021-02-05T00:00:00+00:00

Disable ASLR when debugging with LLDB on macOS

TL;DR

If you want to disable ASLR in LLDB, set the following option:

(lldb) settings set target.disable-aslr false

The long version

This is one of these weird, short posts that the vast majority would not even think of doing since, in most cases, it actually hinders the debugging process of your binary, but here goes. As you might be aware, for a few weeks now, I’m on a mission to write an LLD replacement in Zig, the ZigLd or zld with a first-class support for Mach-O linking and special focus on cross-compilation support (compile anywhere, run on Apple Silicon, or the like). I’ve made some cool progress recently where I am almost able to link a basic (which is not that basic if you take a closer look at) Zig hello-world-style binary. However, when executed I’m greeted with a very informative (not!) segmentation fault message:

> zig build-obj hello.zig
> zld hello.o compiler_rt.o -o hello
> ./hello
[1]    63559 segmentation fault  ./hello

You’d think, just stick it in a debugger to see where you’ve messed up address relocations, right? Well…

> lldb ./hello
(lldb) r
Process 63694 launched: '/Users/kubkon/dev/zld/examples/hello' (arm64)
info: Hello, World!

Process 63694 exited with status = 0 (0x00000000)

Hmm, it works! Wait, what? Why? So it took me a while to realise that when the executable is run via LLDB, LLDB by default turns off Address Space Layout Randomisation (ASLR)! Normally, this is a great thing since your binary is loaded into the memory as if it wasn’t a position-independent executable (PIE). In other words, the address ranges for your local symbols you see when dissecting the disassembled machine code via otool or MachOView.app is preserved. Great! Well, not in our case it turns out. In our case, there is some addressing problem with our binary which only surfaces when the ASLR is on (which BTW is on by default on macOS, and especially on Apple Silicon based Macs). In our case, we’d want to replicate the segfault. Well, fret not, it turns out this can be done with the LLDB, and here’s how. Let’s redo our debugging session but this time we’ll tell LLDB not to disable the ASLR:

> lldb ./hello
(lldb) settings set target.disable-aslr false
(lldb) r
Process 63723 launched: '/Users/kubkon/dev/zld/examples/hello' (arm64)
Process 63723 stopped
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=2, address=0x16b497ff0)
    frame #0: 0x0000000104170270 hello`std.debug.panicExtra + 44
hello`std.debug.panicExtra:
->  0x104170270 +44><: str    x0, [sp, #0x20]
    0x104170274 <+48>: add    x8, sp, #0x54             ; =0x54
    0x104170278 <+52>: str    x8, [sp, #0x18]
    0x10417027c +56><: str    x1, [sp, #0x10]
Target 0: (hello) stopped.

Yay! ASLR has not been disabled and we’ve successfully reproduced the segfault conditions, which means we can hopefully work out what went wrong.

Thread-local storage, LLVM, and Apple Silicon: What can go wrong?

2021-01-21T00:00:00+00:00

Thread-local storage, LLVM, and Apple Silicon: What can go wrong?

The release of Apple Silicon in the form of the M1 chip has definitely stirred things up a lot! M1 boasts some scary performance gains over Intel’s even the most powerful i9 chips, while at the same time using much less energy. But, while the hardware is something to look up to, the decision to switch to the ARM architecture and instruction set, every day unravels more and more problems in existing software. In this short post, I’ll describe the root cause for the thread-local storage (TLS) miscompilation (or should I say, mislinking since this is a linker problem) by the LLVM toolchain on macOS. For context, have a look at this outstanding issue in the Zig compiler: ziglang/zig#7527. I should point out here that Zig’s stage1 (i.e., not-self-hosted) compiler uses LLVM for lowering the intermediate representation into the actual ISA. The self-hosted compiler in Debug mode will by default use our own, in-house incremental linker, whereas in Release mode the plan is to use LLVM to leverage years of really clever micro-optimisations that went into it. (I promise I will write an update about the progress with the incremental Mach-O linker soon-ish, and I will make sure to include all the juicy details that I learnt/discovered along the way. In the meantime, please join me at my FOSDEM21 session where I give an overview of the current state-of-the-art.)

DISCLAIMER: Everything written here is the result of me reverse engineering different Mach-O binaries that I crafted using different tools to try and pinpoint differences between programs that work and those that don’t. For this reason, I don’t expect everything I describe here to always be technically accurate. Unfortunately, the docs on the modern Mach-O are very scarce and so naturally, when writing the incremental linker in Zig, I’ve learnt most of the tricks and built up my understanding of the entire file format by reverse engineering different examples.

TL;DR take everything you read here with a hefty pinch of salt!

The scenario

Take this simple Zig source as an example:

threadlocal var x: usize = 0;

pub fn main() void {
    x = 1;
}

In this blog post, we will investigate how this source maps to the final Mach-O binary, and how a TLS variable x is managed. I should point out here, that currently this example will lead to a segfault on Apple Silicon for the reasons we’ll discuss below.

How is TLS laid out in the final Mach-O binary?

Before we dig deeper into the root cause of TLS miscompilation by the LLD, we should first of all get the basics out of the way. If there is TLS within a Mach-O binary, the flags field within the Mach header must contain the flag:

pub const MH_HAS_TLV_DESCRIPTORS: u32 = 0x800000;

Next, the __DATA segment (or more generally, the read-write segment) will contain two additional sections: __thread_vars and __thread_bss. If you look closely, the former is the threaded equivalent of the __data section, while the latter of the __bss section.

`__thread_vars` section

This section is identified by the flag:

pub const S_THREAD_LOCAL_VARIABLES: u32 = 0x13;

The section itself (or its actual contents) should be 8-byte aligned, and is usually padded with zeros. Here, when the binary is loaded by the dynamic loader dyld, the dyld_stub_binder will write the address of the tlv_get_address symbol, which the compiled machine code responsible for initialising the thread-local variable x, will branch to, but more on this later.

For the dyld_stub_binder to know that it needs to populate the address of tlv_get_address function in the space provided within the __thread_vars section, we use the Dynamic Loader Info (LC_DYLD_INFO_ONLY load command) and in particular, the Binding Info section. Within the Binding Info, we guide the dyld_stub_binder to the right segment and offset within that segment which should equal a cell within the __thread_vars section. I am not going to dig into the actual opcodes used within the Binding Info section to drive the dyld in this blog post. Instead, I’d like to point everyone to a very good resource by Jonathan Levin here.

`__thread_bss` section

This section is also 8-byte aligned and contains zerofill for thread-local storage, much like the __bss section for the static global variables. The section is identified by the flag:

pub const S_THREAD_LOCAL_ZEROFILL: u32 = 0x12;

Also, I should point out that this section doesn’t have data representation within the actual Mach-O file, i.e., its file offset points to the beginning of the binary. However, it is required to point to an unoccupied space within the virtual memory.

How do we initialise a thread-local variable in the compiled machine code?

This is actually more straightforward than you may think. Within the symbol that makes use of the TLS, we call tlv_get_address to initialise the variable. How do we do that? Remember that when the binary was loaded, dyld_stub_binder was run, and fetched the address of tlv_get_address and saved it within the __thread_vars section. Therefore, all that we need to do is to fetch that address and branch to the actual symbol via this address. For the sake of example, assume that we are currently at an address 0x100003CE0, and suppose that the __thread_vars section is at a virtual address of 0x100049180, and that the dyld_stub_binder will store the address of tlv_get_address symbol in the first cell of the section, so at the section’s start address 0x100049180. Then, all we need to do is first of all, locate the memory page where __thread_vars is residing at, to then narrow down to the actual cell, to finally load the address stored in that cell in some register we will branch from. How do we do this with ARM64 ISA?

adrp x0, 70      ; locate the page where `__thread_vars` resides
add  x0, x0, 384 ; locate the start of `__thread_vars` section i.e., `0x100049180`
ldr  x8, [x0, 0] ; dereference the contents of `0x100049180` cell and store it in `x8` register
blr  x8          ; branch with link to address within `x8` register

Just to recap, if the address of the first instruction above is 0x100003CE0, then adrp will load the address of 0x100003CE0 + 70 * 0x1000 := 0x100049000. Note that the result is truncated to the nearest 4KB page. Next, we add 384 (0x180 in hex) to the result, so 0x100049000 + 0x180 := 0x100049180 which is, surprise, surprise, the start address of the __thread_vars section and the address of the first cell of that section. Perfect! Then, we simply derefence the value stored in 0x100049180, which by this time, will be populated by the dyld_stub_binder and will hold the actual address of the tlv_get_address symbol, and store it in x8 register. Finally, we branch to it.

What does the corresponding snippet look like when generated by the LLD?

Unfortunately, the output generated by the LLD uses ldr instead of add instruction leading to a segfault. Why? First of, this is the output generated by the LLD:

adrp x0, 70      
ldr  x0, [x0, 384] ; Oh no! This is incorrect!
ldr  x8, [x0, 0] 
blr  x8

In the above snippet, the use of ldr as the second instruction will first derefence the contents of memory at an address 0x100049000 and then offset that value by 384, whereas, what should happen is the exact opposite. We should offset the address we want to load from by 384 as we did in the original snippet. This will lead to us branching into whatever garbage was stored in 0x100049000 offset by 384 (or 0x180 in hex), and hence, will most likely lead to a segfault or undefined behaviour.

Why is this particularly important for Zig?

I have to admit, I’m using my M1 MBA to drive the development of the incremental linker, and for some time now, I’ve been struggling by not having any stack traces in case I was hitting an assert or panicking somewhere in the codebase. I didn’t investigate the exact cause until very recently, and it turns out, miscompilation of TLS by LLD is to blame. This is because, in Zig, every thread keeps track of its current panic stage. We achieve this via TLS:

// lib/std/debug.zig, line 242
threadlocal var panic_stage: usize = 0;

Therefore, any panic would have to traverse the corrupted codepath that we examined in the snippet above leading to a segfault without printing any stack trace. Let me bring up an example. Consider this simple Zig snippet:

pub fn main() void {
    @panic("NO!");
}

Compiling and running the resulting binary, will unfortunately currently result in:

❯ ./simple
[1]    31670 segmentation fault  ./simple

As an experiment, I decided to hack a manual fix for this by rewriting the offending instruction from ldr to add, and the result is now as expected:

thread 4981258 panic: NO!
/Users/kubkon/dev/zig/examples/st.zig:2:5: 0x10070feff in main (st)
    @panic("NO!");
    ^
/Users/kubkon/dev/zig/build/lib/zig/std/start.zig:335:22: 0x10071004b in std.start.main (st)
        root.main();
                 ^
???:?:?: 0x181b0cf33 in ??? (???)
Panicked during a panic. Aborting.
[1]    31664 abort      ./st

What’s next for Zig?

Clearly, a post-mortem fixup to the binary is not an option since this would require disassembling the entire program searching for any mention of the TLS. An intermediate solution for the time being will be for us to not use TLS on Apple Silicon until the problem is fixed in the LLD itself, or our in-house linker is able to perform well enough as a drop-in replacement.

Stubbing out WASI manually in Rust

2020-04-28T00:00:00+00:00

Stubbing out WASI manually in Rust

Following a really cool blog post “A beginner’s guide to adding a new WASI syscall in Wasmtime” by @RaduM on adding a new WASI syscall to the WASI spec and then stubbing it out in Wasmtime, I thought I’ll contribute to his effort in making Wasm and WASI simpler for the beginners, and delve a bit deeper into stubbing out WASI imports manually when compiling from Rust using the bare target wasm32-unknown-unknown. I guess I’m hopeful this would hint everyone on how wasm-bindgen is doing it, or what is actually happening under-the-hood when compiling directly to WASI in Rust (i.e., using the target wasm32-wasi).

Uhm, so what is it that you’re actually trying to show here?

Excellent question! So here’s the problem we’re trying to solve here. Given a simple, bare Rust lib looking as follows:

fn say_hello() {
  println!("Hello WASI!")
}

how can we configure and stub it out so that when compiled to wasm32-unknown-unknown it’s actually run as a standard executable WASI module? Incidentally, the solution to this problem is pretty straightforward once you know what the WASI spec requires, and how to link to the runtime-provided syscall module. But, one thing at a time.

The solution, step-by-step

Step 0. Get the right tools!

You’ll need:

Wasmtime
Rust toolchain
wasm32-unknown-unknown target installed:

$ rustup target add wasm32-unknown-unknown

After you’ve gathered all three, go ahead and create a Rust library template with:

$ cargo new --lib hello-wasi

Step 1. Configure Cargo.toml

What we actually want here is for the library to compile into a dynamic system library, or cdylib. To do that, we need to tweak Cargo.toml to have the following [lib] section added:

# Cargo.toml
[package]
name = "hello-wasi"
version = "0.1.0"
authors = ["Jakub Konka <kubkon@jakubkonka.com>"]
edition = "2018"

[lib]
crate-type = ["cdylib"] # <-- this is the bit we need to add

[dependencies]

cdylib will allow us to build a dynamically linked Wasm module. With the right entrypoint specified, this will allows to use it as a standard executable Wasm module that you would get by compiling a binary Rust crate to wasm32-wasi target directly. You can find more on this here.

Step 2. Add the required entrypoint

At the time of writing, for a Wasm module to classify as a WASI module, it needs to export an entrypoint called _start. Cool, this is easy, let’s add that in to our lib.rs source:

#[no_mangle]
pub unsafe extern "C" fn _start() {

}

Note the #[no_mangle] attribute which will basically stop Rust compiler from mangling the function’s name, and extern "C" since we need to generate ABI that our Wasm runtime (in this case, Wasmtime) can link to.

Step 3. Print “Hello WASI!” to screen

OK, so this is one is the trickiest step so far. Since we’re stubbing out the WASI imports/exports ourselves, println! macro won’t work as it has no idea where to redirect the output to. If we compiled to wasm32-wasi that would have been taken care of for us, but in this case, we’ll have to do this manually ourselves.

We’ll need to import two WASI syscalls provided by the runtime when executing our module. These are:

fd_write – allows writing to a WASI file descriptor (in our case, we’ll redirect to stdout)
proc_exit – terminates the process with the specified exit code

They’re both provided as part of the wasi_snapshot_preview1 module which we’ll link to in our module. Oh, and btw, you can explore the WASI spec, and in particular, wasi_snapshot_preview1 module in more detail here.

Step 3.1. Declaring the imports

We’ll need the following set of imports declared and types/structs defined:

type Fd = u32;
type Size = usize;
type Errno = i32;
type Rval = u32;

#[repr(C)]
struct Ciovec {
    buf: *const u8,
    buf_len: Size,
}

#[link(wasm_import_module = "wasi_snapshot_preview1")]
extern "C" {
    fn fd_write(fd: Fd, iovs_ptr: *const Ciovec, iovs_len: Size, nwritten: *mut Size) -> Errno;
    fn proc_exit(rval: Rval);
}

A word of explanation what the heck is going on here. In order to link with the syscalls provided by the runtime, we need to match the function name and its signature. For the latter to work out, we need to be careful to properly specify the sizes (in bytes) of the input types, and so, under-the-hood, the runtime is expecting a function fd_write to have the signature

(i32, i32, i32, i32) -> i32

i32 here means that every type is 4 bytes aligned.

Step 3.2. Using the imports to stub out `say_hello`

OK, with this out of the way, we can finally stub out say_hello function with a working replacement of println!:

unsafe fn say_hello() -> Errno {
   let text = "Hello WASI!";
   let ciovec = Ciovec {
       buf: text.as_ptr(),
       buf_len: text.len(),
   };
   let ciovecs = [ciovec];
   let mut nwritten = 0;
   fd_write(1, ciovecs.as_ptr(), ciovecs.len(), &mut nwritten)
}

And our entrypoint:

#[no_mangle]
pub unsafe extern "C" fn _start() {
  let ret = say_hello();
  proc_exit(ret as Rval)
}

And that’s it! We can compile our module and run using Wasmtime.

Step 4. Compile for `wasm32-unknown-unknown`

$ cargo build --target wasm32-unknown-unknown

If there were no errors, you should find your Wasm module in target/wasm32-unknown-unknown/debug/hello_wasi.wasm.

Step 5. Run using `Wasmtime`

$ wasmtime target/wasm32-unknown-unknown/debug/hello_wasi.wasm
Hello WASI!

To check that we’re actually calling the syscalls from our module, we can enable the syscall trace like so:

$ RUST_LOG=wasi_common=trace wasmtime target/wasm32-unknown-unknown/debug/hello_wasi.wasm
 DEBUG wasi_common::ctx > WasiCtx inserting entry PendingEntry::Thunk(0x7ffee24e04c8)
 DEBUG wasi_common::sys::unix::oshandle > Host fd 0 is a char device
 DEBUG wasi_common::sys::unix::oshandle > Host fd 0 is a char device
 DEBUG wasi_common::ctx                 > WasiCtx inserted at Fd(0)
 DEBUG wasi_common::ctx                 > WasiCtx inserting entry PendingEntry::Thunk(0x7ffee24e04c8)
 DEBUG wasi_common::sys::unix::oshandle > Host fd 1 is a char device
 DEBUG wasi_common::sys::unix::oshandle > Host fd 1 is a char device
 DEBUG wasi_common::ctx                 > WasiCtx inserted at Fd(1)
 DEBUG wasi_common::ctx                 > WasiCtx inserting entry PendingEntry::Thunk(0x7ffee24e04c8)
 DEBUG wasi_common::sys::unix::oshandle > Host fd 2 is a char device
 DEBUG wasi_common::sys::unix::oshandle > Host fd 2 is a char device
 DEBUG wasi_common::ctx                 > WasiCtx inserted at Fd(2)
 DEBUG wasi_common::old::snapshot_0::ctx > WasiCtx inserting (0, Some(PendingEntry::Thunk(0x7ffee24e5db0)))
 DEBUG wasi_common::old::snapshot_0::sys::unix::entry_impl > Host fd 0 is a char device
 DEBUG wasi_common::old::snapshot_0::ctx                   > WasiCtx inserting (1, Some(PendingEntry::Thunk(0x7ffee24e5dc0)))
 DEBUG wasi_common::old::snapshot_0::sys::unix::entry_impl > Host fd 1 is a char device
 DEBUG wasi_common::old::snapshot_0::ctx                   > WasiCtx inserting (2, Some(PendingEntry::Thunk(0x7ffee24e5dd0)))
 DEBUG wasi_common::old::snapshot_0::sys::unix::entry_impl > Host fd 2 is a char device
 TRACE wasi_common::wasi::wasi_snapshot_preview1           > fd_write(fd=Fd(1),iovs=*guest 0xfffd8/1)
Hello WASI! TRACE wasi_common::wasi::wasi_snapshot_preview1           >      | result=(nwritten=11)
 TRACE wasi_common::wasi::wasi_snapshot_preview1           >      | errno=No error occurred. System call completed successfully. (Errno::Success(0))
 TRACE wasi_common::wasi::wasi_snapshot_preview1           > proc_exit(rval=0)

Final words

I hope this has been useful for you. If I have made a mistake anywhere (which I most likely have), please feel free to reach out to me via any means available :-)

Porting projects to WASI: the flite program

2019-04-20T00:00:00+00:00

Porting projects to WASI: the flite program

This post is hopefully the first of many more to come demonstrating my efforts to port as many advanced projects written in C/C++ to WASI as possible. What is WASI? WASI is a system interface to run WebAssembly outside the web. If you would like to learn more, I cannot recommend enough the fantastic blog post by Lin Clark Standardising WASI. You should definitely go and check it out!

We’ll start our journey with a fantastic, open-source text-to-speech program called flite. Without further ado, let’s dig in!

Getting the right tools

Firstly, make sure you’ve got all the right tools for the job. In this case, I mean the WASI SDK which you can get from CraneStation/wasi-sdk. The SDK features clang-8 which is capable of targetting wasm32-unknown-wasi triple plus it also features the WASI sysroot. So head over to the website, and get the SDK.

When porting any C program to WASI, there are a few things to watch out for. Firstly, WASI currently doesn’t support setjmp nor sockets so the source code will need to modified accordingly (see here for missing functionality in WASI). In case of flite, it’s somewhat easier as there are two macros DIE_ON_ERROR and CST_NO_SOCKETS which pretty much handle most it for us.

Secondly, remember to correctly set the paths to the compiler, llvm-ar, llvm-ranlib, etc. That is, don’t trust configure to get it right (since WASI is still experimental, the tools won’t be included as OS packages any time soon). Otherwise, you are almost guaranteed to experience missing symbols during linking or other weird errors.

Finally, if the project is using configure, you’ll need to take extra care at including WASI as a supported host. A brilliant, short-and-simple version suggested by Frank Denis is using grep and sed:

$ grep -q -F -- '-wasi' config.sub || sed -i -e 's/-nacl\*)/-nacl*|-wasi)/' config.sub

Porting flite

I’ve modified the source of flite to make it easier to port it to WASI. The source code can be cloned from kubkon/flite. Furthermore, if you’re interested in the changes I’ve had to introduce to the original project, simply compare the fork against the original in Github.

Thus, let’s clone it

$ git clone https://github.com/kubkon/flite

At this point, you need to have the WASI SDK installed. If you don’t have it yet, head over to CraneStation/wasi-sdk and install it. For the rest of this blog post, I’ll assume that you have the SDK installed in /opt/wasi-sdk.

Let’s spin configure and make then!

$ ./configure --host=wasm32-unknown-wasi CC="/opt/wasi-sdk/bin/clang --sysroot=/opt/wasi-sdk/share/sysroot" AR=/opt/wasi-sdk/bin/llvm-ar RANLIB=/opt/wasi-sdk/bin/llvm-ranlib
$ make

If everything went well, you should now have bin/flite program compiled to WASI, and you should be able to run it in any WASI compatible runtime such as CraneStation/wasmtime:

$ wasmtime bin/flite -- --help
flite: a small simple speech synthesizer
  Carnegie Mellon University, Copyright (c) 1999-2016, all rights reserved
  version: flite-2.2-current Sep 2018 (http://cmuflite.org)
usage: flite TEXT/FILE [WAVEFILE]
  Converts text in TEXTFILE to a waveform in WAVEFILE
  If text contains a space the it is treated as a literal
  textstring and spoken, and not as a file name
  if WAVEFILE is unspecified or "play" the result is
  played on the current systems audio device.  If WAVEFILE
  is "none" the waveform is discarded (good for benchmarking)
  Other options must appear before these options
  --version   Output flite version number
  --help      Output usage string
  -o WAVEFILE Explicitly set output filename
  -f TEXTFILE Explicitly set input filename
  -t TEXT     Explicitly set input textstring
  -p PHONES   Explicitly set input textstring and synthesize as phones
  --set F=V   Set feature (guesses type)
  -s F=V      Set feature (guesses type)
  --seti F=V  Set int feature
  --setf F=V  Set float feature
  --sets F=V  Set string feature
  -ssml       Read input text/file in ssml mode
  -b          Benchmark mode
  -l          Loop endlessly
  -voice NAME Use voice NAME (NAME can be pathname/url to flitevox file)
  -voicedir NAME Directory containing (clunit) voice data
  -lv         List voices available
  -add_lex FILENAME add lex addenda from FILENAME
  -pw         Print words
  -ps         Print segments
  -psdur      Print segments and their durations (end-time)
  -pr RelName Print relation RelName
  -voicedump FILENAME Dump selected (cg) voice to FILENAME
  -v          Verbose mode
$ wasmtime --dir=. bin/flite -- test.txt test.wav
$ file test.wav
test.wav: RIFF (little-endian) data, WAVE audio, Microsoft PCM, 16 bit, mono 8000 Hz

I hope this was useful for you! If you have any questions, comments or suggestions, feel free to drop me a line!

The ECB art!

2017-09-12T00:00:00+00:00

The ECB art!

Do you remember from your crypto course why you should never use ECB as a block cipher? Well, here is why:

It looks great though!

The script

In order to compile the code below, you will need to install Crypto++ library, and download header-only bitmap library.

#include <cryptopp/aes.h>
#include <cryptopp/osrng.h>
#include <cryptopp/modes.h>
#include <iostream>
#include <string>
#include <memory>
#include <cmath>
#include "bitmap_image.hpp"

using namespace CryptoPP;

#define AES_KEYLENGTH 32

class FlattenedRGBImage {
  public:
    FlattenedRGBImage(size_t width, size_t height) : _width(width), _height(height) {
      _red = std::unique_ptr<byte[]>(new byte[_width * _height]);
      _blue = std::unique_ptr<byte[]>(new byte[_width * _height]);
      _green = std::unique_ptr<byte[]>(new byte[_width * _height]);
    }

    byte* red() { return _red.get(); }
    byte* blue() { return _blue.get(); }
    byte* green() { return _green.get(); }

  private:
    size_t _width;
    size_t _height;
    std::unique_ptr<byte[]> _red;
    std::unique_ptr<byte[]> _blue;
    std::unique_ptr<byte[]> _green;
};

int quantise(int in, int levels) {
  auto step = (int) 256 / levels;
  auto n = (int) in / step;
  return n * levels;
}

int main(int argc, char* argv[]) {
  if (argc < 2) {
    std::cerr << "You need pass the full path to the .bmp file!" << std::endl;
    return 1;
  }

  std::string path(argv[1]);

  // load bitmap image
  bitmap_image in_image(path);

  if (!in_image) {
    std::cerr << "Error - Failed to open file " << path << std::endl;
    return 1;
  }

  // flatten for encryption
  const auto height = in_image.height();
  const auto width  = in_image.width();
  FlattenedRGBImage plain(width, height);
  auto ptr_red = plain.red();
  auto ptr_blue = plain.blue();
  auto ptr_green = plain.green();

  for (auto y = 0U; y < height; ++y) {
    for (auto x = 0U; x < width; ++x) {
       rgb_t colour;
       in_image.get_pixel(x, y, colour);
       *(ptr_red++) = (byte)quantise((int)colour.red, 4);
       *(ptr_blue++) = (byte)quantise((int)colour.blue, 4);
       *(ptr_green++) = (byte)quantise((int)colour.green, 4);
    }
  }

  // generate random key
  AutoSeededRandomPool rng;
  SecByteBlock key(0x0, AES_KEYLENGTH);
  rng.GenerateBlock(key, key.size());

  FlattenedRGBImage cipher(width, height);

  try {
    // Encrypt using ECB mode
    ECB_Mode<AES>::Encryption enc(key, key.size());
    enc.ProcessData(cipher.red(), plain.red(), height * width);
    enc.ProcessData(cipher.blue(), plain.blue(), height * width);
    enc.ProcessData(cipher.green(), plain.green(), height * width);

    // Write to bitmap image
    ptr_red = cipher.red();
    ptr_blue = cipher.blue();
    ptr_green = cipher.green();
    bitmap_image out_image(width, height);
    for (auto y = 0U; y < height; ++y) {
      for (auto x = 0U; x < width; ++x) {
         rgb_t colour;
         colour.red = (char)*(ptr_red++);
         colour.blue = (char)*(ptr_blue++);
         colour.green = (char)*(ptr_green++);
         out_image.set_pixel(x, y, colour);
      }
    }
    out_image.save_image("ecb_encrypted.bmp");
  }
  catch(const CryptoPP::Exception& e) {
    std::cerr << e.what() << std::endl;
    return 1;
  }

  return 0;
}

C++ short stories: type traits, concepts, and type constraints

2017-09-02T00:00:00+00:00

C++ short stories: type traits, concepts, and type constraints

Suppose you find yourself in a hypothetical situation where you are designing a complex algorithm, for example, an algorithm that requires multiple depth-first searches through a complex tree-like structure. In situations like these, while I do appreciate the power of the debugger, if each node in the tree is a compound data structure (a sequence of some sort, say), it becomes somewhat cumbersome to ‘step into’ a couple hundred times when all you want is to get an idea what the structure looks like, or the path you took, etc. This is especially true when recursive calls/structures are used. Then, the good-old way of printing to screen can be of great help. Not a problem, just need to make sure that every data structure we work with is convertible to string. But what if our algorithm is inherently generic (think templates); i.e., does not care what the primitive data type stored in the tree is until runtime? How do we communicate to the end-developer using our algorithm to ensure that the data types he applies it to are convertible to string? We could test the primitive type at runtime and throw an exception if need be. But runtime errors are difficult to work with and inherently dangerous. We should take a cue from Rust developers and do as much error handling at the compile-time as possible. Furthermore, is it absolutely necessary to require that those data types are convertible to string? After all, it was a requirement only so that we were able to debug the algorithm while developing it, so perhaps the end-developer should not be burdened with it especially if they do not require convertibility to string. So, the puzzle is as follows: given a class template, how do we constrain the template argument to be convertible to string? Do we need to constrain the template argument for the entire definition of the class, or can we limit the constraints to methods that actually require the convertibility to string? Thus, in what follows, we will try to design a class template such that it comprises of a method to_string(), and we will see how we can require for the template argument to be convertible to std::string for the entire class definition, and limited to the to_string() method only.

Idea 1 – using type traits

The first approach relies on the use of type traits and, as such, requires a compiler conformant with C++14 standard. In particular, it makes use of the std::enable_if_t construct which is a helper construct for std::enable_if defined as follows:

template<bool B, class T = void>
struct enable_if;

And as per the C++ reference, if B is true, then the struct std::enable_if has a public member std::enable_if::type equal to T. Thus, we can use std::enable_if to enable class if, at compile-time, the template argument satisfies some template constraint encoded as the predicate B. This suggests the first solution to our problem to look as follows:

template<typename T, typename Enable = void>
class A {};

template<typename T>
class A<T, typename std::enable_if_t<std::is_convertible<T, std::string>::value>> {
public:
  std::string to_string() const {
    return "Class A<>";
  }
};

In the snippet above, we have defined a generic class template A with the second template argument defaulting to void. This class will act as a fallback should our constraint be unmet. The second definition of A acts as a specialisation of the class for all types which are convertible to std::string. This is guaranteed by the call:

std::enable_if_t<std::is_convertible<T, std::string>::value>

Here, std::is_convertible evaluates whether T can be implicitly cast to std::string, and returns true if this is the case. This then becomes the predicate B for std::enable_if, and hence the struct std::enable_if has a member std::enable_if::type. Therefore, if T is convertible to std::string, then the second specialisation of A will be matched by the compiler. Otherwise, if T is not convertible to std::string, then std::enable_if::type is undefined, and hence, the first definition of A will be matched. This approach can be somewhat improved by moving the std::enable_if_t construct into the signature of the method requiring std::string convertibility; namely, the to_string() method. This is precisely what is demonstrated in the snippet below:

template<typename T>
class B {
public:
  template<typename ToString = T>
  typename std::enable_if_t<std::is_convertible<ToString, std::string>::value, std::string>
  to_string() const {
    return "Class B<>";
  }
};

Firstly, note that we have created a new template argument ToString that defaults to T. This is to ensure that only B::to_string() is affected by the constraint, and not the entire class B (try it with the class-wide template argument T and see what happens). We have also changed the call to std::enable_if_t a little bit:

std::enable_if_t<std::is_convertible<ToString, std::string>::value, std::string>

This call now implies that if ToString is convertible to std::string, then std::enable_if::type will equate to std::string, and hence, the entire thing will evaluate to the correct return type for the B::to_string() method; namely, to std::string. Ultimately, if T is not convertible to std::string then rather than not being able to create an object of type B<T>, we simply forbid the caller from using the B<T>::to_string() method.

Bringing it all together, we end up with the short example program:

#include <iostream>
#include <string>
#include <type_traits>

template<typename T, typename Enable = void>
class A {};

template<typename T>
class A<T, typename std::enable_if_t<std::is_convertible<T, std::string>::value>> {
public:
  std::string to_string() const {
    return "Class A<>";
  }
};

template<typename T>
class B {
public:
  template<typename ToString = T>
  typename std::enable_if_t<std::is_convertible<ToString, std::string>::value, std::string>
  to_string() const {
    return "Class B<>";
  }
};

struct Size_t {
  Size_t(size_t v) : value(v) {}

  operator std::string() const {
    return std::to_string(value);
  }

  size_t value;
};

int main() {
  // we can create objects of type A and B from a type not convertible to std::string
  // but we are not allowed to call to_string() method both calls below will error
  // out at compile time
  A<size_t> a1;                             // OK
  B<size_t> b1;                             // OK
  std::cout << a1.to_string() << std::endl; // ERROR!
  std::cout << b1.to_string() << std::endl; // ERROR!

  // however, no problem with Size_t type that is a simple encapsulation of size_t with
  // implemented implicit std::string cast operator
  A<Size_t> a2;                             // OK
  B<Size_t> b2;                             // OK
  std::cout << a2.to_string() << std::endl; // OK
  std::cout << b2.to_string() << std::endl; // OK

  return 0;
}

As discussed above, while A<size_t> and B<size_t> both compile fine, we are not allowed to call A<size_t>::to_string() and B<size_t>::to_string() since size_t is not convertible to std::string. However, everything works OK for type Size_t which we defined as a dummy convenience struct which encapsulates size_t and defines the conversion to std::string (via the operator std::string() const).

Idea 2 – using concepts and type constraints

While Idea 1 is functional and will get the job done, it could be simplified substantially with the use of concepts. Concepts are an upcoming feature of C++20 standard. They allow us to define certain constraints on the type itself which can then be caught at compile-time (essentially, std::enable_if on steroids and more). Since concepts are a work-in-progress, if you want to try out the code described below (and I strongly encourage you to), make sure you have the latest GCC and compile the snippets presented below with -fconcepts flag.

Firstly, we need to define an appropriate concept which will constrain a type T to be castable to std::string. This can be done as follows:

template<typename T>
concept bool CastableToString = requires(T a) {
  { a } -> std::string;
};

Note that each concept is a predicate (hence, bool in the signature). In the body of the concept, we only require that, given an instance of type T, it is possible to convert it (cast it) to std::string. With the concept specified, we can use it to constrain the type for either the entire class or for a subset of methods that would require it. This can be done with the requires keyword. And so, in the former case, this can be done as follows:

template<typename T> requires CastableToString<T>
class C {
public:
  std::string to_string() const {
    return "Class C<>";
  }
};

Here, when instantiating class C<T>, T has to satisfy our concept CastableToString. In order to constrain only some methods within the class, we could do so as follows:

template<typename T>
class D {
public:
  std::string to_string() const requires CastableToString<T> {
    return "Class D<>";
  }
};

In this case, even if T is not convertible to std::string, D<T> can be instantiated just fine; however, D<T>::to_string() will be undefined, and calling it will yield a compile-time error.

Bringing it all together, we end up with the short example program:

#include <iostream>
#include <string>
#include <type_traits>

template<typename T>
concept bool CastableToString = requires(T a) {
  { a } -> std::string;
};

template<typename T> requires CastableToString<T>
class C {
public:
  std::string to_string() const {
    return "Class C<>";
  }
};

template<typename T>
class D {
public:
  std::string to_string() const requires CastableToString<T> {
    return "Class D<>";
  }
};

struct Size_t {
  Size_t(size_t v) : value(v) {}

  operator std::string() const {
    return std::to_string(value);
  }

  size_t value;
};

int main() {
  // errors out since size_t violates the constraint CastableToString
  // but D does not require the same constraint; only D::to_string() call does
  C<size_t> c1;                             // ERROR!
  D<size_t> d1;                             // OK
  std::cout << d1.to_string() << std::endl; // ERROR!

  C<Size_t> c2;                             // OK
  D<Size_t> d2;                             // OK
  std::cout << c2.to_string() << std::endl; // OK
  std::cout << d2.to_string() << std::endl; // OK

  return 0;
}

As outlined above, since size_t is not convertible to std::string, C<size_t> yields an immediate compile-time error. D<size_t> compiles fine; however, we are not allowed to call D<size_t>::to_string(). Finally, everything works OK for type Size_t which we defined as a dummy convenience struct which encapsulates size_t and defines the conversion to std::string (via the operator std::string() const).

Conway's Game of Life

2015-03-15T00:00:00+00:00

Conway’s Game of Life

Since it is the first post I’ve written in a longer while, I’ll keep it short and simple. Lately, in my spare time, I’ve been exploring C# and .NET framework. To this end, I’ve had a look through two great books: C# in a Nutshell by P. Drayton, B. Albahari and T. Neward, and Concurrency in C# Cookbook by S. Cleary. In order to further enhance my understanding of the concepts therein presented, I’ve decided to implement something simple and yet engaging. For that wee project of mine, I’ve decided to go for the Conway’s Game of Life.

The rules of the game

The game was invented in 1970 by a British mathematician, John Conway. In the simplest terms, the game constitutes an imitation (simulation) of life. Within the game, the world is represented by a 2D grid, while each element of the grid, referred to as cell, represents living being. Yes, you guessed correctly, each cell is either alive or dead.

The evolution of life is subject to certain grand rules; that is, whether the cell survives or dies from one stage of the simulation to the next one is determined by the following conditions. The cell survives if 1) it is alive and it has exactly 2 alive neighbours (cells within distance of 1 from the current cell); or 2) it has exactly 3 alive neighbours, regardless of whether it currently is alive or dead. In all other cases, sadly, the cell dies.

Example

The beauty of the game comes from the fact that, depending on the initial state of the world (the distribution of live vs dead cells), the simulation may lead to a totally different (and sometimes very surprising) end state.

To illustrate the beauty of the game, I’ve enclosed a GIF of one particular instance of the game:

If I’ve managed to spark your interest and, like me, you find the simulation weirdly mesmerising, feel free to explore the beauty of the Conway’s Game of Life yourself using my C# app. The source code is available here.

Example of streaming data from the database using Persistent and Conduit libraries in Haskell

2014-01-23T00:00:00+00:00

Example of streaming data from the database using Persistent and Conduit libraries in Haskell

At work, I have to deal with increasingly larger datasets. And by large, I mean a PostgreSQL table with over 7 millions rows of data. If I really wanted to, I could probably be able to load it all into RAM…maybe, I haven’t even tried as it just sounds silly, and especially given the fact that I have quite a few more than one monstrous table like that to deal with. At the same time, I would also like to do as much processing in parallel as possible (you know, put those multi-core CPUs and/or GPUs into action). In short, big data problem.

This problem has made me look into Haskell yet again. After all, Haskell is meant to make it much easier for the programmer to parallelise their algorithms. This book by Simon Marlow is well worth a look when it comes to parallel programming in Haskell, by the way.

There are quite a few packages that make it easy (or easier) to interface with a database in Haskell. Personally, I’ve decided to give Persistent a shot. Persistent is the default database interface used by Yesod Web Framework, works with the majority of popular databases (SQLite, PostgreSQL, MySQL, etc.), and is quite well documented.

In this post, I will provide a (hopefully informative) example of how to use Persistent together with Conduit library. Conduit is a Haskell library that allows streaming data (not necessarily from a database) in constant memory. In other words, with Conduit, it is possible to pull into memory just enough data for processing and not more, thus minimising unnecessary memory use.

The problem

Suppose we are given a database with one massive table of integers, and our objective is to print the records as pairwise tuples to screen. That is, if xs :: [Int] is the list of all the records, we want to zip it with its tail and print it to the screen:

Prelude> mapM_ (putStrLn . show) $ zip xs $ tail xs

Create test SQLite database

Firstly, let’s define the structure of the table using Persistent. Persistent relies on Template Haskell to specify the structure of the database tables. For example, the following code (the source code for all the snippets in this post can be found on GitHub):

-- specs.hs
{-# LANGUAGE EmptyDataDecls    #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE GADTs             #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE QuasiQuotes       #-}
{-# LANGUAGE TemplateHaskell   #-}
{-# LANGUAGE TypeFamilies      #-}

module Specs where

import qualified Database.Persist.TH as TH

TH.share [TH.mkPersist TH.sqlSettings, TH.mkMigrate "migrateAll"] [TH.persistLowerCase|
MyRecord
  value Int
  deriving Show
|]

generates the following SQL:

CREATE TABLE "my_record" ("id" INTEGER PRIMARY KEY, "value" INTEGER NOT NULL);

Having specified the table structure, let’s go ahead and create an SQLite database with a single table with 100000 integers starting from 1 and incrementing by 1:

-- createDB.hs
{-# LANGUAGE OverloadedStrings #-}

import Database.Persist (insertMany)
import Database.Persist.Sqlite (runSqlite,runMigration)
import qualified Specs as S

main :: IO ()
main = runSqlite "test.sqlite" $ do

  -- Create test DB and populate with values
  let n = 100000
  runMigration S.migrateAll
  insertMany $ map S.MyRecord [1..n]

  return ()

When compiled and executed, this code creates a dummy SQLite database “test.sqlite” with a table consisting of 100000 rows of integers.

Print all records as pairs using traditional list approach

In order to query the database for all the records, Persistent provides a helper function, selectList :: [Filter val] -> [SelectOpt val] -> m [Key val]. This function takes two lists containing search criteria as arguments: a list of SQL filter constraints (WHERE), and a list of other SELECT options (e.g., ASC or LIMIT). The result is a list enclosed in a monad m of all the records matching the search criteria. Therefore, in order for us to get all the records from our dummy database, it suffices to call the function with two empty lists:

records <- selectList [] []

All that remains to be done, is to extract the integer values from the list of records into a new list, zip the resultant list with its tail, and print the result to screen. That is,

let values = map (myRecordValue . entityVal) records
mapM_ (liftIO . putStrLn . show) $ zip values $ tail values

You’ll notice three odd functions used in here: myRecordValue, entityVal and liftIO. The first one, myRecordValue is inferred by Template Haskell and represents the constructor of MyRecord data type. If you recall the definition of the SQL table used in this example:

TH.share [TH.mkPersist TH.sqlSettings, TH.mkMigrate "migrateAll"] [TH.persistLowerCase|
MyRecord
  value Int
  deriving Show
|]

Persistent, behind the scenes, translates this template into Haskell data type of the following structure:

data MyRecord = MyRecord
  { myRecordValue :: !Int
  }
  deriving (Show, Read, Eq)

Thus, we can use myRecordValue to access the underlying Int value from MyRecord data type.

The entityVal function, on the other hand, is a constructor for Entity data type, and can be used to access the record data type that’s stored within it. That is, Persistent wraps up every record received from the database in Entity data type. Therefore, in our example, as selectList returns a list [Entity MyRecord] (wrapped in a monad), it is necessary to first recover MyRecord from the Entity, and that can be accomplished via the use of entityVal function.

Lastly, liftIO function is used to lift a computation from the IO monad. In our example, we use it to lift (or complete) an internal IO-type computation used by Persistent, and return back to the standard IO monad.

The full code looks like this:

-- selectList.hs
{-# LANGUAGE OverloadedStrings #-}

import Control.Monad (mapM_)
import Control.Monad.IO.Class  (liftIO)
import Database.Persist (selectList,entityVal)
import Database.Persist.Sqlite (runSqlite)
import Specs (myRecordValue)

main :: IO ()
main = runSqlite "test.sqlite" $ do

    -- Select all records from DB and unpack into a list of Ints
    records <- selectList [] []
    let values = map (myRecordValue . entityVal) records
    mapM_ (liftIO . putStrLn . show) $ zip values $ tail values

If we now run the compiled version of the code with the statistics flag -s; that is,

$ ./selectList +RTS -s

the following output will be generated:

   1,242,454,680 bytes allocated in the heap
      71,087,496 bytes copied during GC
      15,219,032 bytes maximum residency (7 sample(s))
       1,563,864 bytes maximum slop
              35 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0      2379 colls,     0 par    0.05s    0.05s     0.0000s    0.0008s
  Gen  1         7 colls,     0 par    0.03s    0.04s     0.0059s    0.0181s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.97s  (  1.24s elapsed)
  GC      time    0.08s  (  0.10s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    1.05s  (  1.34s elapsed)

  %GC     time       7.6%  (7.1% elapsed)

  Alloc rate    1,277,517,423 bytes per MUT second

  Productivity  92.4% of total user, 72.7% of total elapsed

Notice the total memory used of 35 MB.

Print all records as pairs using Conduit approach

Now, let’s do the same, but this time, let’s use the Conduit library to stream the data from the database. In order to stream the data, we need to have a Conduit Source. Persistent provides a helper function that is much like selectList but returns a Source rather than a list. It’s not difficult to guess its name which is selectSource. The function has a signature selectSource :: [Filter val] -> [SelectOpt val] -> Source m (Entity val).

Having acquired the Source of our data, we now need to extract the underlying Int from the input data type Entity MyRecord, combine the extracted integers into pairwise tuples, convert the tuples into their String representation, and finally, print them to screen.

In order to extract Ints from the upstream, we can use Data.Conduit.List.map helper function. This function applies the supplied function to the data provided by the upstream. In our case, selectSource returns data of type Entity MyRecord, therefore, we can use a composition of myRecordValue with entityVal to extract the underlying Int from input. That is,

import qualified Data.Conduit.List as CL

entityToValue :: Monad m => Conduit (Entity MyRecord) m Int
entityToValue = CL.map (myRecordValue . entityVal)

Note that the type of entityToValue function is Conduit (Entity MyRecord) m Int. This can be read as: for every input of type Entity MyRecord I, the Conduit, will do some processing and pass an Int forward to the downstream (for other Conduits or a Sink to receive).

The next thing to do is to combine the Ints into pairs, and cast them into String. In order to achieve this, we can use the await, yield and leftover helper functions provided by Data.Conduit package. The idea is thus to await for two consecutive Ints, wrap them in a tuple, convert it to String, yield it back to the downstream, and pass the second of the Ints (a leftover) back to the upstream. In code, that would go like this:

showPairs :: Monad m => Conduit Int m String
showPairs = do
  mi1 <- await -- get the next value from the input stream
  mi2 <- await
  case (mi1, mi2) of
    (Just i1, Just i2) -> do
      yield $ show (i1, i2) -- pass tuple of Ints converted
                            -- to String downstream
      leftover i2           -- pass the second component of
                            -- the tuple back to itself (to
                            -- the upstream)
      showPairs
    _ -> return ()

Finally, all that remains is to print each tuple to screen, and terminate the flow of the data stream. To accomplish this, we’ll use the Data.Conduit.List.mapM_ which applies a monadic action to all values in the stream. It’s not difficult to guess that the action we are going to apply is printing the tuples to screen. That is,

printString :: (Monad m, MonadIO m) => Sink String m ()
printString = CL.mapM_ (liftIO . putStrLn)

Having all the components defined, we can link them into the following data stream which will result in the desired outcome:

selectSource [] [] $$ entityToValue =$ showPairs =$ printString

The full code looks like this:

-- selectSource.hs
{-# LANGUAGE OverloadedStrings #-}

import Control.Monad.IO.Class  (MonadIO,liftIO)
import Data.Conduit (await,yield,leftover,Conduit,Sink,($$),(=$))
import qualified Data.Conduit.List as CL
import Database.Persist (selectSource,entityVal,Entity)
import Database.Persist.Sqlite (runSqlite)
import Specs (MyRecord(..),myRecordValue)

-- |Unpacks incoming values from upstream from MyRecord to Int
entityToValue :: Monad m => Conduit (Entity MyRecord) m Int
entityToValue = CL.map (myRecordValue . entityVal)

-- |Converts pairwise tuples of Ints into String
showPairs :: Monad m => Conduit Int m String
showPairs = do
  mi1 <- await -- get the next value from the input stream
  mi2 <- await
  case (mi1, mi2) of
    (Just i1, Just i2) -> do
      yield $ show (i1, i2) -- pass tuple of Ints converted
                            -- to String downstream
      leftover i2           -- pass the second component of
                            -- the tuple back to itself (to
                            -- the upstream)
      showPairs
    _ -> return ()

-- |Prints input String
printString :: (Monad m, MonadIO m) => Sink String m ()
printString = CL.mapM_ (liftIO . putStrLn)

main :: IO ()
main = runSqlite "test.sqlite" $ do

  -- Select all records from DB and return them as a Source
  selectSource [] [] $$ entityToValue =$ showPairs =$ printString

If we now run the compiled version of the code with a statistics flag -s; that is,

$ ./selectSource +RTS -s

the following output will be generated:

   3,152,019,584 bytes allocated in the heap
       7,109,552 bytes copied during GC
          68,928 bytes maximum residency (2 sample(s))
          21,120 bytes maximum slop
               1 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0      6041 colls,     0 par    0.03s    0.04s     0.0000s    0.0000s
  Gen  1         2 colls,     0 par    0.00s    0.00s     0.0003s    0.0006s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    2.70s  (  3.18s elapsed)
  GC      time    0.03s  (  0.04s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    2.74s  (  3.22s elapsed)

  %GC     time       1.3%  (1.2% elapsed)

  Alloc rate    1,167,680,199 bytes per MUT second

  Productivity  98.7% of total user, 83.8% of total elapsed

Note that the total memory used equals only 1 MB compared to 35 MB when using the list approach! I don’t know about you, but I’m impressed.

Unconstrained nonlinear programming: Nelder-Mead Simplex algorithm

2013-10-22T00:00:00+00:00

Unconstrained nonlinear programming: Nelder-Mead Simplex algorithm

In the previous two posts, I have described the basics of penalty methods (exterior and interior), and how they can be used to solve constrained nonlinear optimisation problems using unconstrained nonlinear programming algorithms. In this post, I will take a step back, and describe the simplest (yet quite powerful) heuristic algorithm for unconstrained nonlinear programming: the Nelder-Mead Simplex algorithm.

The algorithm

The algorithm tries to minimise the objective function, \(f\), by successively re-evaluating the so-called simplices in the direction where the function decreases in value [1], [2]. A simplex is a generalisation of the triangle to multi-dimensional spaces. In a 2-dimensional space, the simplex is still just a triangle. Therefore, to keep things simple, the entire post will assume minimisation in 2 dimensions.

Suppose we are given simplex with the following vertices \(\{x_0, x_1, x_2\}\), \(x_i\in\mathbb{R}^2\) for each \(i\). Let’s re-order the vertices by the increasing function value; that is, let \(\{x^L, x^M, x^H\}\) such that \(f(x^L) \le f(x^M) \le f(x^H)\). The main idea of the algorithm is to replace the point which maximises the function in the current simplex with a new, better point, \(x'\) say, in the sense that it minimises the function with respect to \(x^H\); that is, \(f(x') < f(x^H)\). This can be achieved in one of the 4 possible ways: by reflection, expansion, contraction, and shrinking.

Reflection

In the reflection step, the point for which the function achieves the highest value, \(x^H\), is potentially replaced with a point that is reflected against the edge of the simplex connecting the remaining points. That is, if you imagine the edge connecting the points \(x^L\) — \(x^M\) being a mirror, we reflect \(x^H\) against that mirror:

Or mathematically

\[\begin{equation} x^R = \bar{x} + \alpha(\bar{x} - x^H) \end{equation}\]

where \(\alpha>0\) (and usually \(\alpha=1\)), and \(\bar{x} = \frac{x^L + x^M}{2}\) is the centroid (or the centre of gravity) of the points \(x^L\) and \(x^M\). The reasoning behind the reflection step is that since the points \(x^L\) and \(x^M\) minimise the function for the current simplex, by reflecting \(x^H\) against these two points, we should potentially move in the direction that further minimises the function; that is, \(f(x^R) \le f(x^H)\).

Expansion

If it so happens that the reflected point is better than the current minimum, i.e., if \(f(x^R) < f(x^L)\), then perhaps by moving even further in that direction, we’ll further minimise the function. This step is referred to as the expansion step:

Or mathematically

\[\begin{equation} x^E = \bar{x} + \gamma(x^R - \bar{x}) \end{equation}\]

where \(\gamma > 1\), and typically \(\gamma = 2\).

Contraction

If, however, the reflected point, \(x^R\) does not offer a significant improvement over \(x^H\), i.e., if \(f(x^H) \ge f(x^R) > f(x^M)\), then we need to test for points between either the reflected point and the centroid, or the highest value point and the centroid. If \(f(x^R) < f(x^H)\), then we choose the first pair; otherwise, we settle on the latter one. Having selected the appropriate pair of points, we contract the distance between the two. Suppose for the moment that we have selected \(x^R\) to contract against. Then, the contraction step looks as follows:

Or mathematically

\[\begin{equation} x^C = \bar{x} + \beta(x^R - \bar{x}) \end{equation}\]

where \(0 < \beta < 1\), and typically \(\beta = \frac{1}{2}\).

Shrinking

If the contraction step fails, that is, if \(f(x^C) \ge f(x^H)\), then, as the last resort, we shrink the simplex by half towards the current minimum, \(x^L\):

The algorithm revisited

Taking it all together, the algorithm proceeds as follows. For any given simplex \(\{x^L, x^M, x^H\}\), compute the reflected point, \(x^R\).

If \(f(x^R) < f(x^L)\), then compute the expanded point, \(x^E\). If \(f(x^E) < f(x^R)\), then assign \(x^H := x^E\), and start over. Otherwise, assign \(x^H := x^R\), and start over.

If \(f(x^L) \le f(x^R) \le f(x^M)\), then assign \(x^H := x^R\), and start over. Otherwise, select \(x' = \arg\min(f(x^H), f(x^R))\). Compute the contracted point, \(x^C\). If \(f(x^C) \le f(x')\), then assign \(x^H := x^C\), and start over. Otherwise, shrink the simplex, and start over.

Examples

Having described the theory, let’s look at some examples. The full source code for the examples shown in this post can be downloaded from GitHub.

Firstly, consider minimising the function

\[\begin{equation} f(x_1, x_2) = x_1^2 + x_2^2 - 3x_1 - x_1x_2 + 3. \end{equation}\]

This function has one global minimum at a point \(x = (2, 1)\). The path of the simplex algorithm is overlaid on the contour plot of the objective function:

And a more interactive version:

Notice that the algorithm after a series of reflection, expansion, contraction, and shrinking steps converges to the minimum.

Multiple minima

It becomes quite problematic for the algorithm if the function under consideration has more than one minimum. Then, the algorithm is not guaranteed to converge to the global one. Consider an example

\[\begin{equation} f(x_1, x_2) = x_1^4 + x_2^4 - \frac{5}{4} x_1^2 + \frac{1}{4}. \end{equation}\]

This function has two minima at points \((0.8, 0)\) and \((-0.8, 0)\). Depending on the choice of the initial simplex, the algorithm might converge to either of the two. Therefore, it is important to correctly pick the initial simplex for a given problem.

For example, starting off with the simplex \(\{(-2,0), (-1,2), (-1.75,2)\}\) produces the following result:

While, starting off at \(\{(2,0), (1,2), (1.75,2)\}\), produces:

No minimum

Finally, if the function does not have a minimum, then the algorithm will inevitably fail to converge (which is as expected). However, if you are attempting to implement your own version of the Simplex algorithm, it is important to remember about this eventuality. After all, in real world, it is sometimes almost impossible to determine whether the problem we are trying to optimise has an optimum value.

Consider an example

\[\begin{equation} f(x_1, x_2) = x_1x_2. \end{equation}\]

The function does not have a minimum (nor a maximum). However, it does have a stationary point at \((0,0)\). If we run the simplex algorithm in the search of a minimum, it will never converge. The following contour plot shows the search path for 8 iterations of the algorithm:

And a more interactive version:

References

[1] Avriel, M., “Nonlinear Programming: Analysis and Methods”, chapter 9, Dover Publications, 2003.

[2] Mathews, J. H., and Fink, K. K., “Numerical Methods Using Matlab, 4th edition”, chapter 8, Prentice-Hall Inc., 2004. Excerpt on Nelder-Mead Simplex algorithm available here.

Constrained nonlinear programming: interior penalty methods

2013-10-19T00:00:00+00:00

Constrained nonlinear programming: interior penalty methods

In the previous post, I have briefly outlined how to solve a constrained nonlinear problem using unconstrained nonlinear programming; in particular, using exterior penalty methods. In this post, I’ll present an alternative approach called interior penalty methods (or barrier function methods).

The problem revisited

Let’s revisit the optimisation problem introduced in the previous post. Suppose you want to minimise the function \(f(x) = x^2\) subject to the constraint \(x \geq 1\). The function attains its minimum at \(x = 1\); however, if considered without the constraint, the minimum (clearly) shifts to \(x = 0\). Recall that the unconstrained optimisation method will find the latter as the solution to the problem since we can’t enforce the feasible region on the unconstrained optimisation algorithm.

Interior penalty function

However, we can trick the algorithm into converging on the desired solution using the so-called interior penalty function [1]. The function’s aim is to penalise the unconstrained optimisation method if it tries to leave (cross the boundary of) the feasible region of the problem. Hence, with this approach, we are required to start searching for the optimal solution of the problem within the feasible region. This is in direct opposition to the exterior penalty methods where we start outside the feasible region and slowly converge on a minimum from the outside.

Let’s create an interior penalty function for our example:

\[\begin{equation} F(x,\rho^k) = x^2 - \rho^k\ln(x-1) \end{equation}\]

The first derivative of \(F\) with respect to \(x\) looks as follows:

\[\begin{equation} \frac{\partial F}{\partial x} = 2x - \rho^k\cdot\frac{1}{x - 1} \end{equation}\]

Equating the derivative to zero, and noting that by definition \(x > 1\) (since we start in the feasible region of the problem), yields the solution to the modified optimisation problem:

\[\begin{equation} x = \frac{1}{2} + \frac{\sqrt{1 + 2\rho^k}}{2} \end{equation}\]

Note that as \(\rho^k \rightarrow 0\) with \(k\rightarrow\infty\), the solution converges to \(x = 1\). That is, the solution to the original optimisation problem!

This method can readily be converted into a numerical algorithm that uses an unconstrained optimisation method. Pick a number \(\rho\) such that \(0 < \rho < 1\). Starting from \(k = 1\), minimise \(F(x, \rho^k)\) using any unconstrained optimisation method. Finally, feed in the results of the minimisation as the new starting point to the minimisation method, and increment \(k\).

Cython & Python implementation

To finish off this blog post, an exemplary implementation of the interior penalty method in Cython and Python. Like in the previous post, numerical minimisation is performed using algorithms provided by the GNU Scientific Library.

The objective function to be minimised writted in Cython:

from cython_gsl cimport *

from libc.math cimport pow, log

cdef double min_f(double x, void * params) nogil:
    # Extract params
    cdef double * ps = <double *> params
    # Compute interior penalty function
    cdef double barrier = (-1) * ps[0] * log(x - 1)
    # Compute original function value
    cdef double f = pow(x, 2)

    return f + barrier

And C extension to Python that uses Brent’s minimisation method provided by the GNU Scientific Library:

def solve(initial_point, penalty):
    """
    Minimises Cython function min_f using Brent's method, and
    returns the result of the minimisation.

    Arguments:
    initial_point -- initial guess at the minimum
    penalty -- the penalty parameter (rho)
    """
    cdef int status
    # Initialise counter
    cdef int iter = 0
    cdef int max_iter = 100

    # Initial region of convergence
    # NOTE that the region of convergence encompasses
    # all x's greated than 1
    cdef double a = 1.0 + 1e-12
    cdef double b = 100.0

    # Initialise the minimisation problem
    cdef gsl_function F
    cdef double * params = [penalty]
    F.function = &min_f
    F.params = params

    # Initialise Brent's method
    cdef gsl_min_fminimizer_type * T
    cdef gsl_min_fminimizer * s
    T = gsl_min_fminimizer_brent
    s = gsl_min_fminimizer_alloc(T)
    gsl_min_fminimizer_set(s, &F, initial_point, a, b)

    status = GSL_CONTINUE

    # Minimise...
    while (status == GSL_CONTINUE and iter < max_iter):
        iter = iter + 1
        status = gsl_min_fminimizer_iterate(s)

        minimum = gsl_min_fminimizer_x_minimum(s)
        a = gsl_min_fminimizer_x_lower(s)
        b = gsl_min_fminimizer_x_upper(s)

        status = gsl_min_test_interval(a, b, 0.001, 0.0)

        if status == GSL_SUCCESS:
            break

    return minimum

Assuming the Cython module containing min_f and solve functions is called interior, the following Python script demonstrates it in action:

from interior import solve

# Initial minimum
minimum = 3

# Penalty parameter
penalties = map(lambda x: x/10, range(9, 0, -1))

for penalty in penalties:
    # Minimise the modified problem, and feed in the result as
    # the new starting point
    minimum = solve(minimum, penalty)
    
    print("penalty=%.1f, minimum=%f" % (penalty, minimum))

With the following output generated:

penalty=0.9, minimum=1.336660
penalty=0.8, minimum=1.306226
penalty=0.7, minimum=1.274597
penalty=0.6, minimum=1.241620
penalty=0.5, minimum=1.207107
penalty=0.4, minimum=1.170820
penalty=0.3, minimum=1.132455
penalty=0.2, minimum=1.091611
penalty=0.1, minimum=1.047722

Clearly, as the penalty, \(\rho^k\), decreases to zero, the solution to the modified optimisation problem approaches the solution to the original problem.

References

[1] Avriel, M., “Nonlinear Programming: Analysis and Methods”, chapter 12, Dover Publications, 2003.

Constrained nonlinear programming: exterior penalty methods

2013-10-13T00:00:00+00:00

Constrained nonlinear programming: exterior penalty methods

Recently, my PhD research has taken me on an exciting journey through the world of nonlinear programming. Nonlinear programming is a subfield of mathematical optimisation that deals with optimisation problems which are nonlinear in nature. For example, minimising a quadratic function of one real variable is considered as nonlinear since the quadratic function itself is nonlinear.

In this blog post, I thought I’m going to share my thoughts and recent discoveries when it comes to solving constrained nonlinear problems using unconstrained programming. In particular, I want to demonstrate a class of methods falling under the name exterior penalty methods.

The problem

Suppose you want to minimise the function \(f(x) = x^2\) subject to the constraint \(x \geq 1\). The function attains its minimum at \(x = 1\); however, if considered without the constraint, the minimum (clearly) shifts to \(x = 0\). Since we can’t enforce the feasible region (in our case, this would correspond to all \(x\)s larger or equal than \(1\)) on the unconstrained optimisation method, it is necessary to modify the problem in a way that tricks the unconstrained method into finding the minimum within the feasible region.

Exterior penalty function

This can be achieved using the so-called exterior penalty function [1]. The function’s aim is to penalise the unconstrained optimisation method if it converges on a minimum that is outside the feasible region of the problem. Applied to our example, the exterior penalty function modifies the minimisation problem like so:

\[\begin{equation} F(x,\rho^k) = x^2 + \frac{1}{\rho^k}\left[\min(0, x-1)\right]^2 \end{equation}\]

Here, \(\rho^k > 0\) for all \(k\in\mathbb{N}_+\) quantifies the penalty. Note that if \(x\) lies inside the feasible region, then the problem reduces to the original problem; that is, \(F(x, \rho^k) = f(x) = x^2\). It can be shown that as the sequence \((\rho^k), k\in\mathbb{N}_+\) approaches \(0\), the solution to the modified problem approaches the solution to the original problem but from the outside of the feasible region.

This method can readily be converted into a numerical algorithm that uses an unconstrained optimisation method. Pick a number \(\rho\) such that \(0 < \rho < 1\). Starting from \(k = 1\), minimise \(F(x, \rho^k)\) using any unconstrained optimisation method. Feed in the results of the minimisation as the new starting point to the minimisation method, and increment \(k\).

Cython & Python implementation

To finish off this blog post, I’ve included an exemplary implementation of the exterior penalty method for the presented problem in Cython and Python. Why Cython? With Cython, it is very easy to hook into the C-based GNU Scientific Library (which in turn provides the most popular unconstrained optimisation algorithms), and then execute the code from the Python level as a C extension to said language.

The objective function to be minimised writted in Cython:

from cython_gsl cimport *

from libc.math cimport pow

cdef double min_f(double x, void * params) nogil:
    # Extract params
    cdef double * ps = <double *> params
    # Compute penalty
    cdef double penalty = pow(ps[0], ps[1])

    cdef double f = pow(x, 2)

    # If outside feasible region, add the penalty
    if 0 > (x - 1):
        f = f + 1 / penalty * pow(x - 1, 2)

    return f

And C extension to Python that uses Brent’s minimisation method provided by the GNU Scientific Library:

def solve(initial_point, penalty, k):
    """
    Minimises Cython function min_f using Brent's method, and
    returns the result of the minimisation.

    Arguments:
    initial_point -- initial guess at the minimum
    penalty -- the penalty parameter (rho)
    k -- nonnegative natural number
    """
    cdef int status
    # Initialise counter
    cdef int iter = 0
    cdef int max_iter = 100

    # Initial region of convergence
    cdef double a = -100.0
    cdef double b = 100.0

    # Initialise the minimisation problem
    cdef gsl_function F
    cdef double * params = [penalty, k]
    F.function = &min_f
    F.params = params

    # Initialise Brent's method
    cdef gsl_min_fminimizer_type * T
    cdef gsl_min_fminimizer * s
    T = gsl_min_fminimizer_brent
    s = gsl_min_fminimizer_alloc(T)
    gsl_min_fminimizer_set(s, &F, initial_point, a, b)

    status = GSL_CONTINUE

    # Minimise...
    while (status == GSL_CONTINUE and iter < max_iter):
        iter = iter + 1
        status = gsl_min_fminimizer_iterate(s)

        minimum = gsl_min_fminimizer_x_minimum(s)
        a = gsl_min_fminimizer_x_lower(s)
        b = gsl_min_fminimizer_x_upper(s)

        status = gsl_min_test_interval(a, b, 0.001, 0.0)

        if status == GSL_SUCCESS:
            break

    return minimum

Assuming the Cython module containing min_f and solve functions is called exterior, the following Python script demonstrates it in action:

from exterior import solve

# Rho (penalty)
rho = 0.1
# Initial minimum
minimum = 0.0

for k in range(1, 15):
    # Minimise the modified problem, and feed in the result as
    # the new starting point
    minimum = solve(minimum, rho, k)
    
    print("k={}, minimum={}".format(k, minimum))

With the following output generated:

k=1, minimum=0.9090909090909095
k=2, minimum=0.990099009900986
k=3, minimum=0.9990009990009904
k=4, minimum=0.9999000099989902
k=5, minimum=0.9999900001000043
k=6, minimum=0.999999000000997
k=7, minimum=0.9999999000000095
k=8, minimum=0.9999999899999891
k=9, minimum=1.0000000049011502
k=10, minimum=1.0000000049011502
k=11, minimum=1.0000000049011502
k=12, minimum=1.0000000049011502
k=13, minimum=1.0000000049011502
k=14, minimum=1.0000000049011502

It is clear that the algorithm yields the correct solution to the minimisation problem, and the precision of the solution increases with each iteration until the saturation level is reached (after 9th iteration).

References

[1] Avriel, M., “Nonlinear Programming: Analysis and Methods”, chapter 12, Dover Publications, 2003.

Circumventing PDF conversion issues on Manuscript Central

2012-11-07T00:00:00+00:00

Circumventing PDF conversion issues on Manuscript Central

This is a rather short post that is meant to serve as a reminder of what to do when Manuscript Central (MC) rejects the PDF version of your paper generated with pdflatex.

By default, pdflatex generates PDFs in version 1.5, and to the best of my knowledge, MC accepts PDFs up to version 1.4. If the PDF is of newer version than 1.4, the following error appears in MC:

Our apologies, but the file(s) your_paper.pdf (PDF Copy) failed to convert to the appropriate pdf and/or html file for review. You may either try uploading the file(s) again in 20 minutes, or contact the Support Team for further assistance. The Support contact information can be found by clicking the “Get Help Now” link in the upper right corner of this page. Read More …

The easiest way to circumvent this error is to force pdflatex to generate a specific version of the PDF. Suppose we settle for version 1.4. This can be achieved like so:

\documentclass[onecolumn, draftcls, 11pt]{IEEEtran}

\pdfminorversion=4 % tell pdflatex to generate PDF in version 1.4

% The rest of the document...

Now, MC should gladly accept the PDF version of your paper.

Lazy evaluation in Python

2012-09-02T00:00:00+00:00

Lazy evaluation in Python

Recently, while working on my PhD, I’ve started exploring different algorithms for computing and approximating Nash equilibria in games (I’m using a really well-written book, “Algorithmic Game Theory” by Nisan et al., and in order to enhance my understanding of the theory, I’m concurrently having a go at implementing some of the algorithms therein presented – you can find an up-to-date result of my work here: py-game-theory). While at it, I’ve bumped (again) into a concept of lazy evaluation. This has motivated me into checking how ‘lazy’ the Python family of languages is.

But before going into detail of my little investigation, a quick recap of the concept of lazy evaluation. According to Wikipedia, “(…) lazy evaluation is an evaluation strategy which delays the evaluation of an expression until its value is needed and which also avoids repeated evaluations. (…)” Therefore, in a lazily evaluated language, it is possible to create an array of elements, and evaluate each element only when needed, e.g., when computing a sum of all elements in the array, and not sooner. A good example of a language which features lazy evaluation by default is Haskell.

When investigating Python for lazy evaluation, I’ve tested both Python 2 and Py3k (recalling that Py3k has supposedly made some serious improvement in that area). Indeed, it turns out that Py3k is ‘lazier’ than Python 2 (IMHO, yet another point in favour of using Py3k rather than Python 2).

For example, in Python 2, a built-in range function returns a list:

>>> range(5)
[0, 1, 2, 3, 4]

While in Py3k, the same built-in function returns an iterable range object:

>>> range(5)
range(0, 5)
>>> list(range(5))
[0, 1, 2, 3, 4]

Therefore, each number in a sequence is evaluated on demand rather than immediately. This is a very useful feature when dealing with large sequences of numbers, and can save valuable execution time. Three other widely-used built-ins, map, filter, and zip, exhibit similar difference in behaviour between the two versions of Python (try it out!).

A really cool implication of lazy evaluation, is the possibility of creating an infinite (yes, infinite) sequence of numbers. For example, the following code shows how to implement an infinite sequence of Fibonacci numbers in Py3k¹:

def fibonacci():
  a, b = 0, 1
  while True:
    yield a
    a, b = b, a + b

This construct is called a generator function in Python’s nomenclature, and allows for lazy evaluation of the function’s definition.

In order to print the first 10 Fibonacci numbers using our generator function, fibonacci, we first need to create an iterator (or technically, a generator) object by calling fibonacci. Then iterate over some control sequence of 10 elements (e.g., a range object of 10 numbers), and use next built-in function on the created iterator to return successive Fibonacci numbers:

fibb = fibonacci()
for i in range(10):
  print(next(fibb), end=' ')

Remember though that after evaluating the first 10 numbers using fibb, you won’t be able to reprint them with the same generator object, as it resumes where it left off previously. Therefore, you either need to create a new generator object, or save the elements in a container for later use, such as a list; for example, using a list comprehension:

fibb = fibonacci()
fibb_list = [next(fibb) for i in range(10)]
print(fibb_list)

Happy lazy evaluating!

¹Implementation of Fibonacci generator function based on an example from Dive into Python 3 book by Mark Pilgrim.

Watery electronics (or how my Macbook Air survived drowning in water)

2012-08-23T00:00:00+00:00

Watery electronics (or how my Macbook Air survived drowning in water)

A full week ago, I experienced something you dread reading about on the web, but never expect it to happen to you: I managed to spill an entire glass of water (approx. 250 ml) on my 10-month-old 11” Macbook Air (MBA). As you can imagine, I was not overly pleased with the incident, and the fact that I was 1.5 weeks after tonsillectomy, did not help in the least.

However, as it turns out, my distraught was for naught (at least for now). Last night (hence, 6.5 days after the incident), I decided to finally take my chances, and tried turning my MBA on. To my great surprise, it fired up fine, and thus far, does not exhibit any glitches common to electronics (and laptops in particular) exposed to liquids; for example, short-circuited keyboard, or automatic shutdowns. That, alas, does not exclude the possibility of impending corrosion.

An interesting question that comes to mind is whether it was sheer luck (in tandem with a degree of sound judgement and quick-wittedness) on my part that saved the laptop from certain oblivion, or, perhaps, Apple in its ingenuity included a precautionary measure in the form of a moisture sensor that turns the laptop off when the moisture reaches dangerously high levels. Food for thought!

Anyhow, for completeness, I will outline the steps I have taken to save my MBA from drowning. Perhaps it will be of use if not to someone else then to me should I ever relive this dreadful experience.

The incident and my attempts at saving the MBA

When I accidentally knocked my only just refilled glass of water, the water mainly spread across the keyboard of the MBA; however, I believe that some of it got into the opening separating the bottom chassis from the screen (where the vent is located). The laptop was at the time of the incident turned on, and connected to the external power source with a MagSafe adapter. Hence, all the ingredients for a disaster were there.

The steps I took in an attempt at remedying the situation can be succinctly summarized as follows. First and foremost, I disconnected the laptop from the external power source. (The laptop shut itself down; hence, there was no need for me to manually turn it off. Furthermore, I restrained from pressing the power button, and checking whether it is indeed turned off in order to avoid turning the laptop on again, and risking the possibility of a short-circuit.). I turned the laptop sideways to get rid off as much water as possible from the keyboard area. With a handful of paper towels, I wiped out the rest. Finally, I let it dry for a week. The laptop was left in an open position, resting on a few paper towels upside down.

Jakub Konka

Disable ASLR when debugging with LLDB on macOS

Disable ASLR when debugging with LLDB on macOS

TL;DR

The long version

Thread-local storage, LLVM, and Apple Silicon: What can go wrong?

Thread-local storage, LLVM, and Apple Silicon: What can go wrong?

The scenario

How is TLS laid out in the final Mach-O binary?

__thread_vars section

__thread_bss section

How do we initialise a thread-local variable in the compiled machine code?

What does the corresponding snippet look like when generated by the LLD?

Why is this particularly important for Zig?

What’s next for Zig?

Stubbing out WASI manually in Rust

Stubbing out WASI manually in Rust

Uhm, so what is it that you’re actually trying to show here?

The solution, step-by-step

Step 0. Get the right tools!

Step 1. Configure Cargo.toml

Step 2. Add the required entrypoint

Step 3. Print “Hello WASI!” to screen

Step 3.1. Declaring the imports

Step 3.2. Using the imports to stub out say_hello

Step 4. Compile for wasm32-unknown-unknown

Step 5. Run using Wasmtime

Final words

Porting projects to WASI: the flite program

Porting projects to WASI: the flite program

Getting the right tools

Porting flite

The ECB art!

The ECB art!

The script

C++ short stories: type traits, concepts, and type constraints

C++ short stories: type traits, concepts, and type constraints

Idea 1 – using type traits

Idea 2 – using concepts and type constraints

Conway's Game of Life

Conway’s Game of Life

The rules of the game

Example

Example of streaming data from the database using Persistent and Conduit libraries in Haskell

Example of streaming data from the database using Persistent and Conduit libraries in Haskell

The problem

Create test SQLite database

Print all records as pairs using traditional list approach

Print all records as pairs using Conduit approach

Unconstrained nonlinear programming: Nelder-Mead Simplex algorithm

Unconstrained nonlinear programming: Nelder-Mead Simplex algorithm

The algorithm

Reflection

Expansion

Contraction

Shrinking

The algorithm revisited

Examples

Multiple minima

No minimum

References

Constrained nonlinear programming: interior penalty methods

Constrained nonlinear programming: interior penalty methods

The problem revisited

Interior penalty function

Cython & Python implementation

References

Constrained nonlinear programming: exterior penalty methods

Constrained nonlinear programming: exterior penalty methods

The problem

Exterior penalty function

Cython & Python implementation

References

Circumventing PDF conversion issues on Manuscript Central

Circumventing PDF conversion issues on Manuscript Central

Lazy evaluation in Python

Lazy evaluation in Python

Watery electronics (or how my Macbook Air survived drowning in water)

Watery electronics (or how my Macbook Air survived drowning in water)

The incident and my attempts at saving the MBA

`__thread_vars` section

`__thread_bss` section

Step 3.2. Using the imports to stub out `say_hello`

Step 4. Compile for `wasm32-unknown-unknown`

Step 5. Run using `Wasmtime`