show deaths

1 files changed, 121 insertions, 33 deletions

diff --git a/clusters.c b/clusters.c
index f919db8..ef80ee3 100644
--- a/clusters.c
+++ b/clusters.c
@@ -1,52 +1,140 @@
+#include <time.h>
 #include <stdlib.h>
+#include <cairo.h>
 #include <math.h>
+#include <stdlib.h>
+#include <stdint.h>
 #include <stdio.h>
 
 static double
-gen_exp (double lambda)
+sample_exp (double lambda)
 {
     double u = drand48 ();
 
-    return log (1 - u) / (-lambda);
+    return -log (1 - u) / lambda;
+}
+
+#define FALSE 0
+#define TRUE 1
+
+#define WIDTH 6000
+#define HEIGHT 20
+#define N_YEARS 10
+#define DEATHS_PER_YEAR 10
+
+#define BG_COLOR        1, 1, 1
+#define CLUSTER_WIDTH   10.0
+#define CLUSTER_COLOR   0.0, 0.6, 0.3, 0.4
+#define LINE_WIDTH      1.0
+#define LINE_COLOR      0.0, 0.0, 0.0
+#define DEATH_WIDTH     1.0
+#define DEATH_COLOR     0.0, 0.0, 0.0
+#define DEATH_HEIGHT    (CLUSTER_WIDTH / HEIGHT)
+
+static double *
+generate_deaths (int n_years, double deaths_per_year, int *n_deaths)
+{
+    double w;
+    double t;
+    double *result = malloc (10 * sizeof (double));
+    double *first;
+
+    t = 0.0;
+    *n_deaths = 0;
+
+    w = sample_exp (deaths_per_year);
+    t += w;
+    while (t < n_years)
+    {
+        result[(*n_deaths)++] = t;
+        result = realloc (result, ((*n_deaths) + 10) * sizeof (double));
+
+        w = sample_exp (deaths_per_year);
+        t += w;
+    }
+
+    return result;
 }
 
 int
 main ()
 {
-    int cluster = 1;
+    cairo_surface_t *surface =
+        cairo_image_surface_create (CAIRO_FORMAT_ARGB32, WIDTH, HEIGHT);
     int i;
-    double avg_cluster = 0.0;
-    double avg_x;
-    int n_k = 0;
-    int smaller = 0;
-    int n_clusters = 0;
-
-#define N   187250000.0
-#define LAMBDA 2.0
-#define K 1
+    double t = 0.0;
+    cairo_t *cr = cairo_create (surface);
+    double y = HEIGHT / 2.0;
+    double w;
+    int n_deaths;
+    int cluster;
+    double *deaths;
+    int in_cluster;
+    double cluster_x0;
+
+    srand48 (time (0));
+
+    deaths = generate_deaths (N_YEARS, DEATHS_PER_YEAR, &n_deaths);
     
-    for (i = 0; i < N; ++i)
+    cairo_set_source_rgb (cr, BG_COLOR);
+    cairo_paint (cr);
+    
+    cairo_move_to (cr, 0.0, y);
+
+    for (i = 1; i < n_deaths; ++i)
+    {
+        double d0 = deaths[i - 1];
+        double d1 = deaths[i];
+        double x0 = (d0 / N_YEARS) * WIDTH;
+        double x1 = (d1 / N_YEARS) * WIDTH;
+
+        cairo_set_line_width (cr, LINE_WIDTH);
+        cairo_set_source_rgb (cr, LINE_COLOR);
+        cairo_move_to (cr, x0, y);
+        cairo_line_to (cr, x1, y);
+        cairo_stroke (cr);
+    }
+
+    for (i = 0; i < n_deaths; ++i)
+    {
+        double d = deaths[i];
+        double x = (d / N_YEARS) * WIDTH;
+
+        cairo_set_line_width (cr, DEATH_WIDTH);
+        cairo_move_to (cr, x, (0.5 - DEATH_HEIGHT / 2.0) * HEIGHT);
+        cairo_line_to (cr, x, (0.5 + DEATH_HEIGHT / 2.0) * HEIGHT);
+        cairo_set_source_rgb (cr, DEATH_COLOR);
+        cairo_stroke (cr);
+    }
+
+    in_cluster = 0;
+    for (i = 1; i < n_deaths; ++i)
     {
-	double f;
-	f = gen_exp (LAMBDA);
-	if (f <= (1/(LAMBDA)))
-	{
-	    cluster++;
-	    smaller++;
-	}
-	else
-	{
-	    avg_cluster += cluster + 1;
-	    if (cluster == K)
-		n_k++;
-	    cluster = 0;
-	    n_clusters++;
-	}
-	avg_x += f;
+        double d0 = deaths[i - 1];
+        double d1 = deaths[i];
+        double x1 = (d1 / N_YEARS) * WIDTH;
+
+        if (!in_cluster && d1 - d0 <= 1.0 / DEATHS_PER_YEAR)
+        {
+            cluster_x0 = (d0 / N_YEARS) * WIDTH;
+            in_cluster = 1;
+        }
+        if (in_cluster && d1 - d0 > 1.0 / DEATHS_PER_YEAR)
+        {
+            double x1 = (d0 / N_YEARS) * WIDTH;
+
+            cairo_set_line_cap (cr, CAIRO_LINE_CAP_SQUARE);
+            cairo_move_to (cr, cluster_x0, y);
+            cairo_line_to (cr, x1, y);
+            cairo_set_line_width (cr, CLUSTER_WIDTH);
+            cairo_set_source_rgba (cr, CLUSTER_COLOR);
+            cairo_stroke (cr);
+
+            in_cluster = 0;
+        }
     }
+    
+    cairo_surface_write_to_png (surface, "diagram.png");
 
-    printf ("Probability of %d: %f\n", K, n_k / N);
-    printf ("Average cluster size: %f\n", avg_cluster / ((double)n_clusters));
-    printf ("Average value: %f\n", avg_x / N);
-    printf ("Probability of smaller: %f\n", smaller / N);
+    return 0;
 }