Merge branch 'timeline_unit_test' into 'master'

Timeline unit tests

See merge request !1543

.gitignore

@@ -179,6 +179,7 @@ tests/testHashmap
 tests/testNeutrinoCosmology
 tests/testNeutrinoFermiDirac
 tests/testLog
+tests/testTimeline
 tests/*.png
 tests/*.txt
 
@@ -378,4 +379,4 @@ sympy-plots-for-*.tex/
 black_formatting_env
 
 # vscode
-*.json
+*.json
\ No newline at end of file

src/timeline.h

@@ -56,7 +56,7 @@ typedef int8_t timebin_t;
  * @param bin The time bin of interest.
  */
 __attribute__((const)) static INLINE integertime_t
-get_integer_timestep(timebin_t bin) {
+get_integer_timestep(const timebin_t bin) {
 
   if (bin <= 0) return 0;
   return 1LL << (bin + 1);
@@ -71,9 +71,11 @@ get_integer_timestep(timebin_t bin) {
  * We use a fast (but exact for any non-zero value) logarithm in base 2
  * calculation based on the bit representation of the number:
  * log_2(x) = (number of bits in the type) - (number of leading 0-bits in x) - 1
+ *
+ * @param time_step An integer time-step length.
  */
 __attribute__((const)) static INLINE timebin_t
-get_time_bin(integertime_t time_step) {
+get_time_bin(const integertime_t time_step) {
 
   /* ((int) log_2(time_step)) - 1 */
   return (timebin_t)((8 * sizeof(integertime_t) - 2) -
@@ -86,8 +88,8 @@ get_time_bin(integertime_t time_step) {
  * @param bin The time bin of interest.
  * @param time_base the minimal time-step size of the simulation.
  */
-__attribute__((const)) static INLINE double get_timestep(timebin_t bin,
-                                                         double time_base) {
+__attribute__((const)) static INLINE double get_timestep(
+    const timebin_t bin, const double time_base) {
 
   return get_integer_timestep(bin) * time_base;
 }
@@ -95,12 +97,14 @@ __attribute__((const)) static INLINE double get_timestep(timebin_t bin,
 /**
  * @brief Returns the integer time corresponding to the start of the time-step
  * given by a time-bin.
+ * If the current time is a possible beginning for the given time-bin, return
+ * the current time minus the time-step size.
  *
  * @param ti_current The current time on the integer time line.
  * @param bin The time bin of interest.
  */
 __attribute__((const)) static INLINE integertime_t
-get_integer_time_begin(integertime_t ti_current, timebin_t bin) {
+get_integer_time_begin(const integertime_t ti_current, const timebin_t bin) {
 
   const integertime_t dti = get_integer_timestep(bin);
   if (dti == 0)
@@ -110,20 +114,27 @@ get_integer_time_begin(integertime_t ti_current, timebin_t bin) {
 }
 
 /**
- * @brief Returns the integer time corresponding to the start of the time-step
+ * @brief Returns the integer time corresponding to the end of the time-step
  * given by a time-bin.
+ * If the current time is a possible end for the given time-bin, return the
+ * current time.
  *
  * @param ti_current The current time on the integer time line.
  * @param bin The time bin of interest.
  */
 __attribute__((const)) static INLINE integertime_t
-get_integer_time_end(integertime_t ti_current, timebin_t bin) {
+get_integer_time_end(const integertime_t ti_current, const timebin_t bin) {
 
   const integertime_t dti = get_integer_timestep(bin);
   if (dti == 0)
     return 0;
-  else
-    return dti * ceil((double)ti_current / (double)dti);
+  else {
+    const integertime_t mod = ti_current % dti;
+    if (mod == 0)
+      return ti_current;
+    else
+      return ti_current - mod + dti;
+  }
 }
 
 /**
@@ -132,7 +143,7 @@ get_integer_time_end(integertime_t ti_current, timebin_t bin) {
  * @param time The current point on the time line.
  */
 __attribute__((const)) static INLINE timebin_t
-get_max_active_bin(integertime_t time) {
+get_max_active_bin(const integertime_t time) {
 
   if (time == 0) return num_time_bins;
 
@@ -149,7 +160,7 @@ get_max_active_bin(integertime_t time) {
  * @param ti_old The last synchronisation point on the time line.
  */
 __attribute__((const)) static INLINE timebin_t
-get_min_active_bin(integertime_t ti_current, integertime_t ti_old) {
+get_min_active_bin(const integertime_t ti_current, const integertime_t ti_old) {
 
   const timebin_t min_bin = get_max_active_bin(ti_current - ti_old);
   return min_bin;

src/tools.c

@@ -793,7 +793,9 @@ void engine_single_force(double *dim, long long int pid,
 }
 
 /**
- * Returns a random number (uniformly distributed) in [a,b[
+ * Returns a random number (uniformly distributed) in [a,b)
+ *
+ * This function is *not* thread-safe.
  */
 double random_uniform(double a, double b) {
   return (rand() / (double)RAND_MAX) * (b - a) + a;

tests/Makefile.am

@@ -35,7 +35,7 @@ TESTS = testGreetings testMaths testReading.sh testKernel testKernelLongGrav \
 	testCbrt testCosmology testOutputList testFormat.sh \
 	test27cellsStars.sh test27cellsStarsPerturbed.sh testHydroMPIrules \
         testAtomic testGravitySpeed testNeutrinoCosmology.sh testNeutrinoFermiDirac \
-	testLog testDistance
+	testLog testDistance testTimeline
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGrav \
@@ -49,7 +49,7 @@ check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGra
 		 testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
 		 test27cellsStars_subset testCooling testComovingCooling testFeedback testHashmap \
                  testAtomic testHydroMPIrules testGravitySpeed testNeutrinoCosmology \
-		 testNeutrinoFermiDirac testLog
+		 testNeutrinoFermiDirac testLog testTimeline
 
 # Rebuild tests when SWIFT is updated.
 $(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -162,6 +162,8 @@ testHashmap_SOURCES = testHashmap.c
 
 testLog_SOURCES = testLog.c
 
+testTimeline_SOURCES = testTimeline.c
+
 testHydroMPIrules = testHydroMPIrules.c
 
 # Files necessary for distribution

tests/testTimeline.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+#include "timeline.h"
+#include "tools.h"
+
+/* How many times to repeat randomly drawn tests. */
+#define NREPEAT 1048576
+
+/**
+ * @brief test get_integer_time_end() function.
+ * For each particle time bin, pick a random valid time_end for
+ * the dt given by the particle bin; then set a random ti_current
+ * by subtracting some time interval < dt from the expected time_end
+ * and see whether the recovered time_end matches up.
+ *
+ * @param bin_min lowest bin to start test from
+ * @param bin_max highest bin to run test with
+ * @param max_nr_timesteps_test maximal number of timesteps for a sim
+ */
+void test_get_integer_time_end(timebin_t bin_min, timebin_t bin_max,
+                               integertime_t max_nr_timesteps_test) {
+
+  integertime_t dti, step, max_step, set_time_end, ti_current, displacement;
+  integertime_t time_end_recovered;
+
+  for (timebin_t bin = bin_min; bin < bin_max; bin++) {
+
+    dti = get_integer_timestep(bin);
+
+    /* Check random state in timeline */
+    /* ------------------------------ */
+
+    for (int r = 0; r < NREPEAT; r++) {
+      /* First pick a place to set this time_end on the timeline. */
+
+      /* we can't have more than this many steps of this size */
+      max_step = max_nr_timesteps_test / dti;
+      if (max_step == 0)
+        error("max step == 0? bin %d max_nr_steps %lld", bin,
+              max_nr_timesteps_test);
+
+      /* Set the time_end at any step in between there */
+      step = (integertime_t)(random_uniform(1., max_step));
+      set_time_end = step * dti;
+
+      /* Do some safety checks */
+      if (set_time_end % dti != 0) error("time_end not divisible by dti?");
+      if (set_time_end > max_nr_timesteps_test)
+        error("Time end > max_nr_timesteps?");
+      if (set_time_end < (integertime_t)0) error("Time end < 0?");
+
+      /* Now mimick a "current time" by removing a fraction of dti from
+       * the step, and see whether we recover the correct time_end */
+      displacement = (integertime_t)(random_uniform(0., 1. - 1e-12) * dti);
+      ti_current = set_time_end - displacement;
+
+      /* Another round of safety checks */
+      if (ti_current > set_time_end)
+        error(
+            "current>time_end? current=%lld time_end=%lld dti=%lld "
+            "displacement=%lld bin=%d\n",
+            ti_current, set_time_end, dti, displacement, bin);
+      if (ti_current < 0)
+        error(
+            "current<0? current=%lld time_end=%lld dti=%lld "
+            "displacement=%lld bin=%d\n",
+            ti_current, set_time_end, dti, displacement, bin);
+
+      /* Now the actual check. */
+      time_end_recovered = get_integer_time_end(ti_current, bin);
+
+      if (time_end_recovered != set_time_end) {
+        error(
+            "time_end incorrect: expect=%lld got=%lld diff=%lld; current=%lld "
+            "displacement=%lld, dti=%lld, bin=%d",
+            set_time_end, time_end_recovered, set_time_end - time_end_recovered,
+            ti_current, displacement, dti, bin);
+      }
+    }
+
+    /* Check time_end = 0 */
+    /* ------------------ */
+    set_time_end = 0;
+    ti_current = 0;
+    time_end_recovered = get_integer_time_end(ti_current, bin);
+    if (time_end_recovered != set_time_end) {
+      error(
+          "time_end incorrect: expect=%lld got=%lld diff=%lld; current=%lld "
+          "bin=%d",
+          set_time_end, time_end_recovered, set_time_end - time_end_recovered,
+          ti_current, bin);
+    }
+
+    /* Check time_end = final_step */
+    /* --------------------------- */
+    set_time_end = max_nr_timesteps_test;
+
+    if (dti < NREPEAT) {
+      /* Check all possible time states before the end */
+      for (integertime_t delta = 1; delta < dti; delta++) {
+        ti_current = max_nr_timesteps_test - delta;
+        time_end_recovered = get_integer_time_end(ti_current, bin);
+        if (time_end_recovered != set_time_end) {
+          error(
+              "time_end incorrect: expect=%lld got=%lld diff=%lld; "
+              "current=%lld bin=%d",
+              set_time_end, time_end_recovered,
+              set_time_end - time_end_recovered, ti_current, bin);
+        }
+      }
+    } else {
+      /* Draw random states again */
+      for (int r = 0; r < NREPEAT; r++) {
+        /* Do some safety checks */
+        if (set_time_end % dti != 0) error("time_end not divisible by dti?");
+
+        /* Now mimick a "current time" by removing a fraction of dti from
+         * the step, and see whether we recover the correct time_end */
+        displacement = (integertime_t)(random_uniform(0., 1. - 1e-12) * dti);
+        ti_current = set_time_end - displacement;
+
+        /* Another round of safety checks */
+        if (ti_current > set_time_end)
+          error(
+              "current>time_end? current=%lld time_end=%lld dti=%lld "
+              "displacement=%lld bin=%d\n",
+              ti_current, set_time_end, dti, displacement, bin);
+        if (ti_current < 0)
+          error(
+              "current<0? current=%lld time_end=%lld dti=%lld "
+              "displacement=%lld bin=%d\n",
+              ti_current, set_time_end, dti, displacement, bin);
+
+        /* Now the actual check. */
+        time_end_recovered = get_integer_time_end(ti_current, bin);
+
+        if (time_end_recovered != set_time_end) {
+          error(
+              "time_end incorrect: expect=%lld got=%lld diff=%lld; "
+              "current=%lld "
+              "displacement=%lld, dti=%lld, bin=%d",
+              set_time_end, time_end_recovered,
+              set_time_end - time_end_recovered, ti_current, displacement, dti,
+              bin);
+        }
+      }
+    }
+  }
+
+  printf("Passed %s\n", __func__);
+}
+
+/**
+ * @brief test get_time_bin by converting all possible bins
+ * into timesteps and trying to recover the initial bin
+ * by calling get_time_bin()
+ * @param bin_min lowest bin to start test from
+ * @param bin_max highest bin to run test with
+ *
+ **/
+void test_get_time_bin(timebin_t bin_min, timebin_t bin_max) {
+
+  for (timebin_t bin = bin_min; bin < bin_max; bin++) {
+    integertime_t dti = get_integer_timestep(bin);
+    timebin_t bin_recovered = get_time_bin(dti);
+    if (bin != bin_recovered) {
+      error("Expected bin=%d, got=%d", bin, bin_recovered);
+    }
+  }
+
+  printf("Passed %s\n", __func__);
+}
+
+/**
+ * @brief test get_integer_time_begin() function.
+ * For each particle time bin, pick a random valid time_begin for
+ * the dt given by the particle bin; then set a random ti_current
+ * by adding some time interval < dt to the expected time_begin
+ * and see whether the recovered time_begin matches up.
+ *
+ * @param bin_min lowest bin to start test from
+ * @param bin_max highest bin to run test with
+ * @param max_nr_timesteps_test maximal number of timesteps for a sim
+ */
+void test_get_integer_time_begin(timebin_t bin_min, timebin_t bin_max,
+                                 integertime_t max_nr_timesteps_test) {
+
+  integertime_t dti, step, max_step, set_time_begin, ti_current, displacement;
+  integertime_t time_begin_recovered;
+
+  for (timebin_t bin = bin_min; bin < bin_max; bin++) {
+
+    dti = get_integer_timestep(bin);
+
+    /* Check random state in timeline */
+    /* ------------------------------ */
+
+    for (int r = 0; r < NREPEAT; r++) {
+      /* First pick a place to set this time_end on the timeline. */
+
+      /* we can't have more than this many steps of this size */
+      max_step = max_nr_timesteps_test / dti;
+      if (max_step == 0)
+        error("max step == 0? bin %d max_nr_steps %lld", bin,
+              max_nr_timesteps_test);
+
+      /* Set the time_end at any step in between there */
+      step = (integertime_t)(random_uniform(0., max_step));
+      set_time_begin = step * dti;
+
+      /* Do some safety checks */
+      if (set_time_begin % dti != 0)
+        error("set time_begin not divisible by dti?");
+      if (set_time_begin >= max_nr_timesteps_test)
+        error("Time begin %lld >= max_nr_timesteps %lld?", set_time_begin,
+              max_nr_timesteps);
+      if (set_time_begin < (integertime_t)0) error("Time begin < 0?");
+
+      /* Now mimick a "current time" by removing a fraction of dti from
+       * the step, and see whether we recover the correct time_end */
+      displacement = (integertime_t)(random_uniform(0., 1.) * dti);
+      ti_current = set_time_begin + displacement;
+
+      /* Another round of safety checks */
+      if (ti_current < set_time_begin)
+        printf(
+            "current<time_begin? current=%lld time_end=%lld dti=%lld "
+            "displacement=%lld bin=%d\n",
+            ti_current, set_time_begin, dti, displacement, bin);
+
+      /* Now the actual check. */
+      time_begin_recovered = get_integer_time_begin(ti_current, bin);
+
+      if (ti_current == set_time_begin) {
+        /* If the displacement was zero, then the function shall return
+         * the beginning of the timestep that ends at ti_current */
+        if (time_begin_recovered + dti != set_time_begin)
+          error(
+              "time_begin incorrect: expect=%lld got=%lld diff=%lld; "
+              "current=%lld "
+              "displacement=%lld, dti=%lld",
+              /*expect=*/set_time_begin - dti,
+              /*got=*/time_begin_recovered,
+              /*diff=*/set_time_begin - dti - time_begin_recovered, ti_current,
+              displacement, dti);
+      } else {
+        if (time_begin_recovered != set_time_begin)
+          error(
+              "time_begin incorrect: expect=%lld got=%lld diff=%lld; "
+              "current=%lld displacement=%lld, dti=%lld",
+              set_time_begin, time_begin_recovered,
+              set_time_begin - time_begin_recovered, ti_current, displacement,
+              dti);
+      }
+    }
+
+    /* Check time_begin = 0 */
+    /* ------------------ */
+    set_time_begin = 0;
+    ti_current = 0;
+    time_begin_recovered = get_integer_time_begin(ti_current, bin);
+    if (time_begin_recovered != set_time_begin)
+      error(
+          "time_begin incorrect: expect=%lld got=%lld diff=%lld; "
+          "current=%lld dti=%lld",
+          set_time_begin, time_begin_recovered,
+          set_time_begin - time_begin_recovered, ti_current, dti);
+
+    /* Check time_begin = final_step */
+    /* ----------------------------- */
+    /* This is a tad nonsensial since in this scenario, we're
+     * at an integer time > max_nr_timesteps */
+    set_time_begin = max_nr_timesteps_test;
+
+    if (dti < NREPEAT) {
+      /* Check all possible time states before the end */
+      for (integertime_t delta = 1; delta < dti; delta++) {
+        ti_current = max_nr_timesteps_test + delta;
+        time_begin_recovered = get_integer_time_begin(ti_current, bin);
+        if (time_begin_recovered != set_time_begin)
+          error(
+              "time_begin incorrect: expect=%lld got=%lld diff=%lld; "
+              "current=%lld displacement=%lld, dti=%lld",
+              set_time_begin, time_begin_recovered,
+              set_time_begin - time_begin_recovered, ti_current, displacement,
+              dti);
+      }
+    } else {
+      /* Draw random states again */
+      for (int r = 0; r < NREPEAT; r++) {
+
+        /* Now mimick a "current time" by removing a fraction of dti from
+         * the step, and see whether we recover the correct time_end */
+        displacement = (integertime_t)(random_uniform(0., 1.) * dti);
+        ti_current = set_time_begin + displacement;
+
+        /* Another round of safety checks */
+        if (ti_current < set_time_begin)
+          error(
+              "current>time_begin? current=%lld time_begin=%lld dti=%lld "
+              "displacement=%lld bin=%d\n",
+              ti_current, set_time_begin, dti, displacement, bin);
+        if (ti_current < 0)
+          error(
+              "current<0? current=%lld time_begin=%lld dti=%lld "
+              "displacement=%lld bin=%d\n",
+              ti_current, set_time_begin, dti, displacement, bin);
+
+        /* Now the actual check. */
+        time_begin_recovered = get_integer_time_begin(ti_current, bin);
+
+        if (ti_current == set_time_begin) {
+          /* If the displacement was zero, then the function shall return
+           * the beginning of the timestep that ends at ti_current */
+          if (time_begin_recovered + dti != set_time_begin)
+            error(
+                "time_begin incorrect: expect=%lld got=%lld diff=%lld; "
+                "current=%lld "
+                "displacement=%lld, dti=%lld",
+                /*expect=*/set_time_begin - dti,
+                /*got=*/time_begin_recovered,
+                /*diff=*/set_time_begin - dti - time_begin_recovered,
+                ti_current, displacement, dti);
+        } else {
+          if (time_begin_recovered != set_time_begin)
+            error(
+                "time_begin incorrect: expect=%lld got=%lld diff=%lld; "
+                "current=%lld displacement=%lld, dti=%lld",
+                set_time_begin, time_begin_recovered,
+                set_time_begin - time_begin_recovered, ti_current, displacement,
+                dti);
+        }
+      }
+    }
+  }
+
+  printf("Passed %s\n", __func__);
+}
+
+/**
+ * @brief test get_max_active_bin by randomly choosing current times
+ * and manually checking whether a higher bin could be active
+ *
+ * @param bin_max highest bin to run test with
+ * @param max_nr_timesteps_test maximal number of timesteps for a sim
+ **/
+void test_get_max_active_bin(timebin_t bin_max,
+                             integertime_t max_nr_timesteps_test) {
+
+  integertime_t ti_current, dti;
+  timebin_t max_active_bin, testbin;
+
+  /* Test random ti_currents */
+  /* ----------------------- */
+  for (int r = 0; r < NREPEAT; r++) {
+    /* test time 0 later, so use time >= 1 */
+    ti_current = (integertime_t)(random_uniform(1., max_nr_timesteps_test));
+
+    max_active_bin = get_max_active_bin(ti_current);
+    testbin = 0;
+    dti = get_integer_timestep(max_active_bin);
+    for (timebin_t bin = max_active_bin; bin < bin_max && dti <= ti_current;
+         bin++) {
+      if (ti_current % dti == 0) testbin = bin;
+      /* update dti here for exit condition */
+      dti = get_integer_timestep(bin + 1);
+    }
+    if (testbin > max_active_bin)
+      error("Found higher active bin: time=%lld max_active_bin=%d found=%d",
+            ti_current, max_active_bin, testbin);
+  }
+
+  /* Test first 2^bin_max integertimes */
+  /* --------------------------------- */
+  for (timebin_t bin = 1; bin < bin_max; bin++) {
+    ti_current = get_integer_timestep(bin);
+    max_active_bin = get_max_active_bin(ti_current);
+
+    if (max_active_bin != bin)
+      error("Got wrong max_active_bin: Expect=%d got=%d time=%lld", bin,
+            max_active_bin, ti_current);
+  }
+
+  /* Test final 2^bin_max integertimes */
+  /* --------------------------------- */
+  for (timebin_t bin = 1; bin < bin_max; bin++) {
+    dti = get_integer_timestep(bin);
+    ti_current = max_nr_timesteps_test - dti;
+    if (ti_current == 0)
+      error("Testing bin which only fits in once into entire timeline");
+    max_active_bin = get_max_active_bin(ti_current);
+
+    if (max_active_bin != bin) {
+      error("Got wrong max_active_bin: Expect=%d got=%d time=%lld", bin,
+            max_active_bin, ti_current);
+    }
+  }
+
+  /* Make sure everything is active at time zero */
+  /* ------------------------------------------- */
+  max_active_bin = get_max_active_bin(0);
+  /* Note: this should be num_time_bins, not bin_max! */
+  if (max_active_bin != num_time_bins)
+    error("Didn't get max_active_bin = num_time_bins at t=0");
+
+  printf("Passed %s\n", __func__);
+}
+
+/**
+ * @brief test get_min_active_bin by randomly choosing current times
+ * for a pair of small and big bins
+ *
+ * @param bin_min smallest bin to run test with
+ * @param bin_max highest bin to run test with
+ * @param max_nr_timesteps_test maximal number of timesteps for a sim
+ **/
+void test_get_min_active_bin(timebin_t bin_min, timebin_t bin_max,
+                             integertime_t max_nr_timesteps_test) {
+
+  integertime_t dti_lo, dti_hi;
+  integertime_t step, max_step;
+  integertime_t ti_old, ti_current;
+  timebin_t min_active_bin;
+
+  /* Test random ti_currents */
+  /* ----------------------- */
+
+  for (timebin_t bin_lo = bin_min; bin_lo < bin_max - 1; bin_lo++) {
+    dti_lo = get_integer_timestep(bin_lo);
+
+    for (timebin_t bin_hi = bin_lo + 1; bin_hi < bin_max; bin_hi++) {
+      dti_hi = get_integer_timestep(bin_hi);
+      max_step = (max_nr_timesteps_test / dti_hi);
+
+      if (NREPEAT / bin_max < max_step) {
+        /* Do random draws */
+        for (int r = 0; r < NREPEAT / bin_max; r++) {
+
+          step = (integertime_t)(random_uniform(0., max_step));
+          ti_current = step * dti_hi;
+          ti_old = ti_current - dti_lo;
+
+          if (ti_current % dti_hi != 0)
+            error("Time current not divisible by dti_hi");
+          if (ti_current % dti_lo != 0)
+            error("Time current not divisible by dti_lo");
+          if (ti_current <= ti_old)
+            error("ti_current=%lld <= ti_old=%lld | bins %d %d; ", ti_current,
+                  ti_old, bin_lo, bin_hi);
+
+          min_active_bin = get_min_active_bin(ti_current, ti_old);
+          if (min_active_bin != bin_lo)
+            error("Got wrong min active bin. Expect=%d got=%d", bin_lo,
+                  min_active_bin);
+        }
+      } else {
+        /* systematically check every possibility */
+        for (step = 0; step <= max_step; step++) {
+
+          ti_current = step * dti_hi;
+          ti_old = ti_current - dti_lo;
+
+          if (ti_current % dti_hi != 0)
+            error("Time current not divisible by dti_hi");
+          if (ti_current % dti_lo != 0)
+            error("Time current not divisible by dti_lo");
+          if (ti_current <= ti_old)
+            error("ti_current=%lld <= ti_old=%lld | bins %d %d; ", ti_current,
+                  ti_old, bin_lo, bin_hi);
+
+          min_active_bin = get_min_active_bin(ti_current, ti_old);
+          if (min_active_bin != bin_lo)
+            error("Got wrong min active bin. Expect=%d got=%d", bin_lo,
+                  min_active_bin);
+        }
+      }
+    }
+  }
+
+  printf("Passed %s\n", __func__);
+}
+
+/**
+ * @brief Check the timeline functions.
+ */
+int main(int argc, char* argv[]) {
+
+  /* Initialize CPU frequency, this also starts time. */
+  unsigned long long cpufreq = 0;
+  clocks_set_cpufreq(cpufreq);
+
+  /* Get some randomness going */
+  const int seed = time(NULL);
+  message("Seed = %d", seed);
+  srand(seed);
+
+  /* run test with the limits we use in SWIFT */
+  /* ---------------------------------------- */
+  test_get_time_bin(1, num_time_bins);
+  test_get_integer_time_begin(1, num_time_bins, max_nr_timesteps);
+  test_get_integer_time_end(1, num_time_bins, max_nr_timesteps);
+  test_get_max_active_bin(num_time_bins, max_nr_timesteps);
+  test_get_min_active_bin(1, num_time_bins, max_nr_timesteps);
+
+  /* run test beyond the limits we use in SWIFT */
+  /* ------------------------------------------ */
+  /* Get maximal bin and number of timesteps based on type size */
+  size_t timebin_bits = sizeof(timebin_t) * 8;
+  timebin_t max_num_time_bins = 0;
+  for (size_t i = 0; i < timebin_bits - 1; i++) {
+    max_num_time_bins |= (1 << i);
+  }
+  size_t timestep_bits = sizeof(integertime_t) * 8;
+
+  /* timestep_bits -1 because of how << works,
+   * additional -1 because integertime_t is signed
+   * additional -1 so we can do any timestep at least twice */
+  max_num_time_bins = min((size_t)max_num_time_bins, timestep_bits - 3);
+
+  if (num_time_bins < max_num_time_bins) {
+    /* Use analogous definition as in timeline.h here */
+    integertime_t max_nr_timesteps_test = (1LL << (max_num_time_bins + 1));
+
+    test_get_time_bin(num_time_bins, max_num_time_bins);
+    test_get_integer_time_begin(num_time_bins, max_num_time_bins,
+                                max_nr_timesteps_test);
+    test_get_integer_time_end(num_time_bins, max_num_time_bins,
+                              max_nr_timesteps_test);
+    test_get_max_active_bin(max_num_time_bins, max_nr_timesteps_test);
+    test_get_min_active_bin(1, max_num_time_bins, max_nr_timesteps_test);
+  }
+
+  return 0;
+}