Cimba

A multithreaded discrete event simulation library in C

What is it?

A fast discrete event simulation library written in C and assembly with POSIX pthreads. Simulated processes are implemented as stackful coroutines ("fibers") inside the pthreads.

Implementation status:

• x86-64: Stable, both for Linux and Windows
• Apple Silicon: Planned
• ARM: Planned

Cimba models run 40-50 times faster than SimPy equivalents. The chart below shows the number of simulated events processed per second of wall clock time on a simple M/M/1 queue implemented in SimPy and Cimba. Cimba runs this scenario 45 times faster than SimPy with all CPU cores in use. Cimba runs 25 % faster (20M events/sec) on a single core than SimPy using all 64 cores (16M events/sec).

Why should I use it?

It is fast, powerful, reliable, and free.

• Fast: The speed from multithreaded parallel execution translates to high resolution in your simulation modeling. You can run hundreds of replications and parameter variations in just a few seconds, generating tight confidence intervals in your experiments and a high density of data points along parameter variations.

  In the benchmark shown above, Cimba reduces the run time by 97.8 % compared to the same model in SimPy using all CPU cores. This translates into doing your simulation experiments in seconds instead of minutes, or in minutes instead of hours.

• Powerful: Cimba provides a comprehensive toolkit for discrete event simulation:

  • Processes implemented as asymmetric stackful coroutines. A simulated process can yield and resume control from any level of a function call stack, allowing well-structured coding of arbitrarily large simulation models. This makes it natural to express agentic behavior by conceptually placing oneself "inside" that process and describing what it does. A simulated process can run in an infinite loop or act as a one-shot customer passing through the system, being both an active agent and a passive object as needed.

  • Pre-packaged process interaction mechanisms like resources, resource pools, buffers, object queues, priority queues, and timeouts. Cimba also provides condition variables where your simulated processes can wait for arbitrarily complex conditions to become true – anything you can express as a function returning a binary true or false result.

  • A wide range of fast, high-quality random number generators, both of academically important and more empirically oriented types. Important distributions like normal and exponential are implemented by state-of-the-art ziggurat rejection sampling for speed and accuracy.

  • Integrated logging and data collection features that make it easy to get a model running and understand what is happening inside it, including custom asserts to pinpoint sources of errors.

  • As a C library, Cimba allows easy integration with other libraries and programs. You could call CUDA routines to enhance your simulation models with GPU-powered agentic behavior or drive a fancy graphics interface like a 3D visualization of a manufacturing plant. You could even call the Cimba simulation engine from other programming languages, since the C calling convention is standard and well-documented.

• Reliable: Cimba is well-engineered open source. There is no mystery to the results you get. The code is written with liberal use of assertions to enforce preconditions, invariants, and postconditions in each function. The assertions act as self-enforcing documentation on expected inputs to and outputs from the Cimba functions. About 13 % of all code lines in the Cimba library are asssertions, a very high density. There are unit tests for each module. Running the unit test battery in debug mode (all assertions active) verifies the correct operation in great detail. You can do that by the one-liner meson test -C build from the terminal command line.

• Free: Cimba should fit well into the budget of most research groups.

What can I use Cimba for?

It is a general-purpose discrete event simulation library, in the spirit of a 21st century Simula67 descendant. You can use it to model, e.g.

• computer networks,
• transportation networks,
• operating system task scheduling,
• manufacturing systems and job shops,
• military command and control systems,
• hospital and emergency room patient flows,
• queuing systems like bank tellers and store checkouts,
• urban systems like public transport and garbage collection,
• and quite a few more application domains of similar kinds, where overall system complexity arises from interactions between relatively simple components.

If you look under the hood, you will also find additional reusable internal components. Cimba contains stackful coroutines doing their own thing on thread-safe cactus stacks. There are fast memory pool allocators for generic small objects and hash-heaps combining a binary heap and an open addressing hash map using Fibonacci hashing. Although not part of the public Cimba API, these components can also be used in your model if needed.

What does the code look like?

It is C11/C17. As an illustration, this is the entire code for our multithreaded M/M/1 benchmark mentioned above:

    #include <inttypes.h>
    #include <stdio.h>
    #include <stdint.h>
    
    #include <cimba.h>
    
    #define NUM_OBJECTS 1000000u
    #define ARRIVAL_RATE 0.9
    #define SERVICE_RATE 1.0
    #define NUM_TRIALS 100
    
    CMB_THREAD_LOCAL struct cmi_mempool objectpool = CMI_MEMPOOL_STATIC_INIT(8u, 512u);
    
    struct simulation {
        struct cmb_process *arrival;
        struct cmb_process *service;
        struct cmb_objectqueue *queue;
    };
    
    struct trial {
        double arr_mean;
        double srv_mean;
        uint64_t obj_cnt;
        double sum_wait;
        double avg_wait;
    };
    
    struct context {
        struct simulation *sim;
        struct trial *trl;
    };
    
    void *arrivalfunc(struct cmb_process *me, void *vctx)
    {
        cmb_unused(me);
        const struct context *ctx = vctx;
        struct cmb_objectqueue *qp = ctx->sim->queue;
        const double mean_hld = ctx->trl->arr_mean;
        for (uint64_t ui = 0; ui < NUM_OBJECTS; ui++) {
            const double t_hld = cmb_random_exponential(mean_hld);
            cmb_process_hold(t_hld);
            void *object = cmi_mempool_alloc(&objectpool);
            double *dblp = object;
            *dblp = cmb_time();
            cmb_objectqueue_put(qp, object);
        }
    
        return NULL;
    }
    
    void *servicefunc(struct cmb_process *me, void *vctx)
    {
        cmb_unused(me);
        const struct context *ctx = vctx;
        struct cmb_objectqueue *qp = ctx->sim->queue;
        const double mean_srv = ctx->trl->srv_mean;
        uint64_t *cnt = &(ctx->trl->obj_cnt);
        double *sum = &(ctx->trl->sum_wait);
        while (true) {
            void *object = NULL;
            cmb_objectqueue_get(qp, &object);
            const double *dblp = object;
            const double t_srv = cmb_random_exponential(mean_srv);
            cmb_process_hold(t_srv);
            const double t_sys = cmb_time() - *dblp;
            *sum += t_sys;
            *cnt += 1u;
            cmi_mempool_free(&objectpool, object);
        }
    }
    
    void run_trial(void *vtrl)
    {
        struct trial *trl = vtrl;
    
        cmb_logger_flags_off(CMB_LOGGER_INFO);
        cmb_random_initialize(cmb_random_hwseed());
        cmb_event_queue_initialize(0.0);
        struct context *ctx = malloc(sizeof(*ctx));
        ctx->trl = trl;
        struct simulation *sim = malloc(sizeof(*sim));
        ctx->sim = sim;
    
        sim->queue = cmb_objectqueue_create();
        cmb_objectqueue_initialize(sim->queue, "Queue", CMB_UNLIMITED);
    
        sim->arrival = cmb_process_create();
        cmb_process_initialize(sim->arrival, "Arrival", arrivalfunc, ctx, 0);
        cmb_process_start(sim->arrival);
        sim->service = cmb_process_create();
        cmb_process_initialize(sim->service, "Service", servicefunc, ctx, 0);
        cmb_process_start(sim->service);
    
        cmb_event_queue_execute();
    
        cmb_process_stop(sim->service, NULL);
        cmb_process_terminate(sim->arrival);
        cmb_process_terminate(sim->service);
        cmb_process_destroy(sim->arrival);
        cmb_process_destroy(sim->service);
    
        cmb_objectqueue_destroy(sim->queue);
        cmb_event_queue_terminate();
        free(sim);
        free(ctx);
    }
    
    int main(void)
    {
        struct trial *experiment = calloc(NUM_TRIALS, sizeof(*experiment));
        for (unsigned ui = 0; ui < NUM_TRIALS; ui++) {
            struct trial *trl = &experiment[ui];
            trl->arr_mean = 1.0 / ARRIVAL_RATE;
            trl->srv_mean = 1.0 / SERVICE_RATE;
            trl->obj_cnt = 0u;
            trl->sum_wait = 0.0;
        }
    
        cimba_run_experiment(experiment,
                             NUM_TRIALS,
                             sizeof(*experiment),
                             run_trial);
    
        struct cmb_datasummary summary;
        cmb_datasummary_initialize(&summary);
        for (unsigned ui = 0; ui < NUM_TRIALS; ui++) {
            const double avg_tsys = experiment[ui].sum_wait / (double)(experiment[ui].obj_cnt);
            cmb_datasummary_add(&summary, avg_tsys);
        }
    
        const unsigned un = cmb_datasummary_count(&summary);
        if (un > 1) {
            const double mean_tsys = cmb_datasummary_mean(&summary);
            const double sdev_tsys = cmb_datasummary_stddev(&summary);
            const double serr_tsys = sdev_tsys / sqrt((double)un);