# Statistical simulation

I've written the following program as a submission to the German Federal Computer Science Competition for highschool students. The program makes the following simulation with a given number of days and a given α. In order to simplify my code, I use a shellscript to run multiple simulations in parallel. You can get all of the code from GitHub: (in directory 5).

I know that my implement correctly the fact that the happiness is calculated once a day. I calculate it just when needed.

## Question

(Translated from German, see source)

The CEO of Nash Ltd. wants to see her workers happy. Thus, she likes to provide them with free coffee. [...] After a period of ten days where the CEO's assistant makes coffee, the workers have to make coffee on their own. Every evening, the remaining coffee is poured away.

The carafe can contain up to ten cups of coffee, there are 15 workers. Each worker likes to drink three cups a day, one in the morning (8 to 10 AM), one at noon (12 to 2 PM) and one in the afternoon (3 to 5 PM). I each of these timeslots the workers arrive in a random order and behave like this:

If there is coffee left in the carafe they pour one cup and leave. If the carafe is empty, the behavior depends on whether they are happy or not. If they are, they brew a new carafe and take one cup. If not, the leave in anger without drinking any coffee.

Whether a worker is happy on a certain day is decided right after they get up. It depends on how often they drank coffee in the last ten days and on how often they brew coffee in that time. The exact definition reads:

If a person drank n cups in the last ten days and brew coffee for m times, he or she is happy iff n >= 10 and m/n < α; α is a personal threshold.

Write a program that simulates the worker where α is equal for all of them. How many cups of coffee is a worker drinking on the long-term average every day depending on α?

## stats.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>

#define STAT_DAYS    10 /* days in statistical window */
#define INTRO_DAYS   10 /* Amount of days the assistant refills the can */
#define WORKER_COUNT 15 /* number of worker in the company */
#define CAN_CAPACITY 10 /* Cups per can */
#define SERVINGS      3 /* Coffee servings per worker per day */

/* data for one worker */
typedef struct {
int cups;
int cans;
} worker_t[STAT_DAYS];

/* data for the whole company at a given point */
typedef struct {
double threshold;
int day_number;
int total_cups; /* total number of cups drunk by the workers */
int can_state;
worker_t workers[WORKER_COUNT];
} company_t;

/* xorshift random number generator for random number generation with good
* statistical properties. See http://en.wikipedia.org/wiki/Xorshift */
static uint32_t xor_state[4], xor_now;
static int xor_entropy = 0;

static uint32_t xor128(int bits) {
uint32_t t;

/* use remaining entropy if possible */
if (xor_entropy < bits) {
xor_entropy = 32;
t = xor_state[0] ^ (xor_state[1] << 11);
xor_state[0] = xor_state[1];
xor_state[1] = xor_state[2];
xor_state[2] = xor_state[3];
xor_state[3] ^= (xor_state[3] >> 19) ^ (t ^ ( t >> 8));
xor_now = xor_state[3];
}

xor_entropy -= bits;
t = xor_now & ((1<<bits)-1);
xor_now >>= bits;
return t;
}

/* initialize the random number generator. Return 0 on success. */
static uint32_t init_xor128(void) {
FILE *entropy = fopen("/dev/urandom","rb");

if (entropy == NULL) {
perror("Can't open /dev/random");
return 1;
}

setbuf(entropy,NULL); /* We don't need any buffering */

fclose(entropy);

fprintf(stderr,"Can't initialize rng.\n");
return 1;
}

return 0;
}

/* number of bits used by i. __builtin_clz is an intrisic available on gcc
* and clang. */
static int bits(uint32_t i) {
return 32 - __builtin_clz(i);
}

/* Generate a permuted array of pointers to members of an array of workers.
* Algorithm taken from http://en.wikipedia.org/wiki/Fisher–Yates_shuffle */
static void shuffle_array(int out[WORKER_COUNT]) {
int i,j;

out[0] = 0;

for (i = 1; i < WORKER_COUNT; i++) {
/* get random numbers till you get one in range */
while (j = xor128(bits(i)), j > i);
out[i] = out[j];
out[j] = i;
}
}

/* Is a worker happy? */
static bool happy(worker_t *worker,double threshold) {
int i;
int cup_count = 0, can_count = 0;

/* get amount of cups drunk, cans brewed */
for (i = 0; i < STAT_DAYS; i++) cup_count += (*worker)[i].cups;
for (i = 0; i < STAT_DAYS; i++) can_count += (*worker)[i].cans;

/* I interprete n >= 10 as "n >= sliding window" */
return cup_count >= STAT_DAYS && 1.0L*can_count/cup_count < threshold;
}

/* simulate one round of coffee for all workers */
static void coffee_round(company_t *company) {
int order[WORKER_COUNT];
int i, day_mod = company->day_number % STAT_DAYS;

shuffle_array(order);

for (i = 0; i < WORKER_COUNT; i++) {
int j = order[i];
if (!company->can_state && company->day_number >= INTRO_DAYS) {
if (!happy(&company->workers[j],company->threshold)) continue;

company->can_state = CAN_CAPACITY - 1;
company->workers[j][day_mod].cans++;
}

company->can_state--;
company->workers[j][day_mod].cups++;
company->total_cups++;
}
}

/* simulate one day of coffee */
static void day(company_t *company) {
/* day_mod: index into the statistical tables */
int i, day_mod = company->day_number % STAT_DAYS;

/* zero out statistics entry for today */
for (i = 0;i < WORKER_COUNT; i++) {
company->workers[i][day_mod].cups = 0;
company->workers[i][day_mod].cans = 0;
}

company->can_state = 0; /* empty can */

for (i = 0;i < SERVINGS; i++) coffee_round(company);

company->day_number++; /* A day has passed */
}

/* Main function */
int main(int argc,char **argv) {
company_t *company;
int day_count;

if (argc != 3) {
fprintf(stderr,"Usage: %s days threshold\n",argv[0]);
return EXIT_FAILURE;
}

/* allocate and set to zero */
company = calloc(sizeof *company,1);

if (!sscanf(argv[1],"%d",&day_count) || day_count <= 0) {
fprintf(stderr,"%s is not a valid number of days.\n",argv[1]);
return EXIT_FAILURE;
}

if (!sscanf(argv[2],"%lf",&company->threshold)) {
fprintf(stderr,"%s is not a valid threshold.\n",argv[2]);
return EXIT_FAILURE;
}

if (init_xor128()) return EXIT_FAILURE;

while (company->day_number < day_count) day(company);

/* print output */
printf("%.6lf\t%.6lf\n",
company->threshold,
((double)company->total_cups/day_count)/WORKER_COUNT);

free(company);

return EXIT_SUCCESS;
}


## stats.sh

#!/bin/bash

#shell script for statistics creation. This shell script generates a table of
#statistics for a lot of happiness thresholds and saves them into happiness.csv

#various parameters
STATFILE=stats.tsv
SIM=./stats
JOBS=$(grep -ic ^processor /proc/cpuinfo) #number of processors available #search space START=0.0 STEP=0.0016 END=0.5 DAYS=$((1024*1024))

LC_NUMERIC=C #set numeric locale to C in order to get correct separators
LC_COLLATE=C #set collating sequence to C in order to sort correctly
#calculate results in parallel
seq $START$STEP $END | xargs -n 1 -P$JOBS $SIM$DAYS >$STATFILE sort -o$STATFILE $STATFILE #sort output as it might be unordered  ## Makefile CFLAGS ?= -O3 -Wall XELATEX ?= xelatex GNUPLOT ?= gnuplot PYGMENTIZE ?= pygmentize PYGMENTIZE_STYLE ?= borland PYGMENTIZE_OPTS ?= -O linenos 5.pdf: 5.tex plot.tex stats.tex definitions.tex$(RM) *.aux
$(XELATEX) -halt-on-error$< >/dev/null
$(XELATEX) -halt-on-error$< >/dev/null

definitions.tex:
$(PYGMENTIZE) -f tex -S$(PYGMENTIZE_STYLE) > $@ .c.tex:$(PYGMENTIZE) $(PYGMENTIZE_OPTS) -o$@ $< plot.tex: stats.gp stats.tsv$(GNUPLOT) $< stats.tsv: stats ./stats.sh clean:$(RM) *.pdf *.aux *.log definitions.tex stats.tex plot.tex *.tsv stats *.eps

.PHONY: clean


Overall, this is very clean and I have only a few small improvements.

1) In happy:

/* get amount of cups drunk, cans brewed */
for (i = 0; i < STAT_DAYS; i++) cup_count += (*worker)[i].cups;
for (i = 0; i < STAT_DAYS; i++) can_count += (*worker)[i].cans;


could be a single for loop.

2) In day (by the way, this should probably be simulate_day or something):

for (i = 0;i < WORKER_COUNT; i++) {
company->workers[i][day_mod].cups = 0;
company->workers[i][day_mod].cans = 0;
}


could be

for (i = 0; i < WORKER_COUNT; i++) {
company->workers[i][day_mod] = {
.cups = 0,
.cans = 0;
}
}


3) In main, you calloc and free a variable; why not just allocate it on the stack?

/* Main function */
int main(int argc,char **argv) {

/* allocate and set to zero */
company = calloc(sizeof *company,1);

/* print output */
printf("%.6lf\t%.6lf\n",


describe what the code is doing, not why it's doing it.

5) There are a few small style inconsistencies (especially with regards to spaces after commas and semicolons). Consider using astyle and writing an .astylerc file for yourself.

• Thank you really much for the feedback! I allocated the structure on the heap as it might be bigger than a few hundred bytes. People always suggest that one should allocate on the heap everything that is not really small. – FUZxxl Dec 26 '12 at 10:13
• Oh. I suppose if WORKER_COUNT goes up, that makes sense. – raylu Dec 27 '12 at 13:12