make mliAvl more useable, turn ground_grid into a histogram. No longe…

…r store the bin assignment for each photon.

libs/mli/src/mliAvlDict.c

@@ -1,4 +1,4 @@
-#include "mliAvlTree.h"
+#include "mliAvlDict.h"
 #include "chk.h"
 #include "mli_math.h"
 
@@ -210,3 +210,10 @@ int mliAvlDict_get(struct mliAvlDict* dict, const int64_t key, int64_t* value)
                 return 1;
         }
 }
+
+void mliAvlDict_reset(struct mliAvlDict* dict)
+{
+        dict->tree.root = NULL;
+        dict->back = 0;
+        dict->len = 0;
+}

libs/mli/src/mliAvlDict.h

@@ -24,5 +24,6 @@ int mliAvlDict_set(struct mliAvlDict* dict, const int64_t key, const int64_t val
 int mliAvlDict_pop(struct mliAvlDict* dict, const int64_t key);
 int mliAvlDict_has(struct mliAvlDict* dict, const int64_t key);
 int mliAvlDict_get(struct mliAvlDict* dict, const int64_t key, int64_t* value);
+void mliAvlDict_reset(struct mliAvlDict* dict);
 
 #endif

libs/mli/src/mliAvlTree.c

@@ -1,4 +1,5 @@
 #include "mliAvlTree.h"
+#include <stdio.h>
 
 /* Adopted from:
  *

libs/mli/src/mli_version.h

@@ -6,7 +6,7 @@
 
 #define MLI_VERSION_MAYOR 1
 #define MLI_VERSION_MINOR 9
-#define MLI_VERSION_PATCH 2
+#define MLI_VERSION_PATCH 3
 
 void mli_logo_fprint(FILE *f);
 void mli_authors_and_affiliations_fprint(FILE *f);

libs/mli_corsika/apps/ground_grid.main.c

@@ -3,6 +3,7 @@
 #include "../src/mli_corsika_ray_voxel_overlap.h"
 #include "../src/mli_corsika_CorsikaPhotonBunch.h"
 #include "../src/mli_corsika_EventTape.h"
+#include "../src/mli_corsika_Histogram2d.h"
 #include "../../mli/src/chk.h"
 #include "../../mli/src/mli_math.h"
 #include "../../mli/src/mliIo.h"
@@ -251,6 +252,7 @@ int main(int argc, char *argv[])
 {
         FILE *istream = NULL;
         FILE *ostream = NULL;
+        struct mliIo buff = mliIo_init();
 
         struct mliEventTapeReader arc = mliEventTapeReader_init();
         float runh[273] = {0.0};
@@ -262,7 +264,8 @@ int main(int argc, char *argv[])
         unsigned int num_bins_each_axis = 1000;
         uint64_t numover = (uint64_t)(1.5 * (double)num_bins_each_axis);
 
-        struct mliGroundGrid grid = mliGroundGrid_init();
+        /*struct mliGroundGrid grid = mliGroundGrid_init();*/
+        struct mliCorsikaHistogram2d hist = mliCorsikaHistogram2d_init();
 
         struct mliCorsikaPhotonBunch bunch;
 
@@ -295,12 +298,27 @@ int main(int argc, char *argv[])
         chk(mliEventTapeReader_begin(&arc, istream));
         chk(mliEventTapeReader_read_runh(&arc, runh));
 
-        chk(mliGroundGrid_malloc(&grid, x_bin_edges.size, y_bin_edges.size));
+        /* chk(mliGroundGrid_malloc(&grid, x_bin_edges.size, y_bin_edges.size));*/
+        chk(mliCorsikaHistogram2d_malloc(&hist, 10*1000));
+        chk(mliIo_malloc_capacity(&buff, 10 * 1000));
 
         while (mliEventTapeReader_read_evth(&arc, evth)) {
                 uint64_t bunch_index = 0;
                 uint64_t xx, yy, ii, jj = 0;
-                mliGroundGrid_reset(&grid);
+                /* mliGroundGrid_reset(&grid); */
+
+                if (hist.dict.capacity > 10 * 1000 * 1000) {
+                        chk(mliCorsikaHistogram2d_malloc(&hist, 10*1000));
+                } else {
+                        mliCorsikaHistogram2d_reset(&hist);
+                }
+
+                if (buff.capacity > 1000 * 1000) {
+                        chk(mliIo_malloc_capacity(&buff, 1000 * 1000));
+                } else {
+                        buff.pos = 0;
+                        buff.size = 0;
+                }
 
                 while (mliEventTapeReader_read_cherenkov_bunch(&arc, raw)) {
                         unsigned int oo;
@@ -315,18 +333,41 @@ int main(int argc, char *argv[])
                                 &z_idxs,
                                 &overlaps);
 
+                        /*
                         chk(mliGroundGrid_assign_overlap(
                                 &grid,
                                 &x_idxs,
                                 &y_idxs,
                                 bunch_index,
                                 bunch.weight_photons));
+                        */
+                        for (oo = 0; oo < overlaps.size; oo ++)
+                        {
+                                chk(
+                                        mliCorsikaHistogram2d_assign(
+                                                &hist,
+                                                x_idxs.array[oo],
+                                                y_idxs.array[oo],
+                                                bunch.weight_photons
+                                        )
+                                );
+                        }
 
                         bunch_index += 1;
                 }
-                chk(mliGroundGrid_fprint_jsonl(&grid, ostream));
+                /*chk(mliGroundGrid_fprint_jsonl(&grid, ostream));*/
+                chk(mliCorsikaHistogram2d_dumps(&hist, &buff));
+                buff.pos = 0;
+                chk_fwrite(
+                        buff.cstr,
+                        sizeof(unsigned char),
+                        buff.size,
+                        ostream
+                );
         }
-        mliGroundGrid_free(&grid);
+        /* mliGroundGrid_free(&grid); */
+        mliCorsikaHistogram2d_free(&hist);
+
         chk(mliEventTapeReader_finalize(&arc));
         fclose(istream);
         fclose(ostream);

libs/mli_corsika/src/mli_corsika_Histogram2d.c

@@ -0,0 +1,144 @@
+/* Copyright 2017 Sebastian A. Mueller */
+
+#include "mli_corsika_Histogram2d.h"
+#include <string.h>
+#include "../../mli/src/chk.h"
+
+struct key {
+        int32_t x;
+        int32_t y;
+};
+
+union i4i4_to_i8 {
+        struct key i4i4;
+        int64_t i8;
+};
+
+struct mliCorsikaHistogram2d mliCorsikaHistogram2d_init(void) {
+        struct mliCorsikaHistogram2d hist;
+        hist.dict = mliAvlDict_init();
+        return hist;
+}
+
+double interpret_int64_as_double(int64_t i) {
+        double f;
+        memcpy(&f, &i, sizeof(double));
+        return f;
+}
+
+int64_t interpret_double_as_int64(double d) {
+        int64_t i;
+        memcpy(&i, &d, sizeof(int64_t));
+        return i;
+}
+
+void mliCorsikaHistogram2d_free(struct mliCorsikaHistogram2d *hist)
+{
+        mliAvlDict_free(&hist->dict);
+}
+
+int mliCorsikaHistogram2d_malloc(
+        struct mliCorsikaHistogram2d *hist,
+        const uint64_t capacity)
+{
+        return mliAvlDict_malloc(&hist->dict, capacity);
+}
+
+int mliCorsikaHistogram2d_assign(
+        struct mliCorsikaHistogram2d *hist,
+        const int32_t x,
+        const int32_t y,
+        const double weight)
+{
+        int has;
+        union i4i4_to_i8 key;
+        int64_t ival = 0;
+        key.i4i4.x = x;
+        key.i4i4.y = y;
+
+        has = mliAvlDict_get(&hist->dict, key.i8, &ival);
+        if (has) {
+                double dval = interpret_int64_as_double(ival);
+                dval += weight;
+                ival = interpret_double_as_int64(dval);
+
+        } else {
+                ival = interpret_double_as_int64(weight);
+        }
+        return mliAvlDict_set(
+                &hist->dict,
+                key.i8,
+                ival
+        );
+}
+
+int mliCorsikaHistogram2d_dumps__(const struct mliAvlNode* node, struct mliIo* f)
+{
+        int64_t count = 0;
+        union i4i4_to_i8 key;
+        double dval = 0.0;
+
+        if (node == NULL) {
+                return 1;
+        }
+
+        key.i8 = node->key;
+        dval = interpret_int64_as_double(node->value);
+
+        count = mliIo_write(
+                f, (const void*)(&key.i4i4.x), sizeof(uint32_t), 1
+        );
+        chk_msg(count == 1, "Failed to write x.");
+
+        count = mliIo_write(
+                f, (const void*)(&key.i4i4.y), sizeof(uint32_t), 1
+        );
+        chk_msg(count == 1, "Failed to write y.");
+
+        count = mliIo_write(
+                f, (const void*)(&dval), sizeof(double), 1
+        );
+        chk_msg(count == 1, "Failed to write weight.");
+
+        if (node->avl.left != NULL) {
+                struct mliAvlNode* left = (struct mliAvlNode*)(node->avl.left);
+                chk_msg(mliCorsikaHistogram2d_dumps__(left, f), "1");
+        }
+        if (node->avl.right != NULL) {
+                struct mliAvlNode* right = (struct mliAvlNode*)(node->avl.right);
+                chk_msg(mliCorsikaHistogram2d_dumps__(right, f), "2");
+        }
+
+        return 1;
+chk_error:
+        return 0;
+}
+
+int mliCorsikaHistogram2d_dumps(
+        const struct mliCorsikaHistogram2d *hist,
+        struct mliIo* f)
+{
+        int64_t count = 0;
+
+        count = mliIo_write(
+                f, (const void*)(&hist->dict.len), sizeof(uint64_t), 1
+        );
+        chk_msg(count == 1, "Failed to write dict->len.");
+
+        chk_msg(
+                mliCorsikaHistogram2d_dumps__(
+                        (const struct mliAvlNode*)hist->dict.tree.root,
+                        f
+                ),
+                "woot?"
+        )
+
+        return 1;
+chk_error:
+        return 0;
+}
+
+void mliCorsikaHistogram2d_reset(struct mliCorsikaHistogram2d *hist)
+{
+        mliAvlDict_reset(&hist->dict);
+}

libs/mli_corsika/src/mli_corsika_Histogram2d.h

@@ -0,0 +1,34 @@
+/* Copyright 2018-2020 Sebastian Achim Mueller */
+#ifndef MLI_CORSIKA_HISTOGRAM2D_H_
+#define MLI_CORSIKA_HISTOGRAM2D_H_
+
+#include <stdint.h>
+#include "../../mli/src/mliAvlDict.h"
+#include "../../mli/src/mliIo.h"
+
+
+struct mliCorsikaHistogram2d{
+        struct mliAvlDict dict;
+};
+
+struct mliCorsikaHistogram2d mliCorsikaHistogram2d_init(void);
+
+int mliCorsikaHistogram2d_malloc(
+        struct mliCorsikaHistogram2d *hist,
+        const uint64_t capacity);
+
+void mliCorsikaHistogram2d_free(struct mliCorsikaHistogram2d *hist);
+
+int mliCorsikaHistogram2d_assign(
+        struct mliCorsikaHistogram2d *hist,
+        const int32_t x,
+        const int32_t y,
+        const double weight);
+
+int mliCorsikaHistogram2d_dumps(
+        const struct mliCorsikaHistogram2d *hist,
+        struct mliIo* f);
+
+void mliCorsikaHistogram2d_reset(struct mliCorsikaHistogram2d *hist);
+
+#endif

libs/mli_corsika/src/mli_corsika_Histogram2d.test.c

@@ -0,0 +1,58 @@
+/* Copyright 2021 Sebastian A. Mueller */
+#include "../../mli_testing/src/mli_testing.h"
+
+CASE("mli_corsika_Histogram2d: init,malloc,free")
+{
+        struct mliCorsikaHistogram2d hist = mliCorsikaHistogram2d_init();
+        CHECK(hist.dict.capacity == 0);
+
+        CHECK(mliCorsikaHistogram2d_malloc(&hist, 100));
+        CHECK(hist.dict.capacity == 100);
+
+        mliCorsikaHistogram2d_free(&hist);
+        CHECK(hist.dict.capacity == 0);
+}
+
+
+CASE("mli_corsika_Histogram2d: assign to same bin multiple times.")
+{
+        uint64_t i = 0;
+        uint64_t a = 0;
+        int32_t b = 0;
+        double c = 0;
+        int64_t rc = 0;
+
+        struct mliCorsikaHistogram2d hist = mliCorsikaHistogram2d_init();
+        struct mliIo buff = mliIo_init();
+
+
+        CHECK(mliCorsikaHistogram2d_malloc(&hist, 10));
+        CHECK(mliIo_malloc(&buff));
+
+        for (i = 0; i < 20; i ++) {
+                CHECK(mliCorsikaHistogram2d_assign(&hist, 42, 13, 0.5));
+        }
+
+        CHECK(mliCorsikaHistogram2d_dumps(&hist, &buff));
+        mliCorsikaHistogram2d_free(&hist);
+
+        buff.pos = 0;
+
+        rc = mliIo_read(&buff, (const void *)(&a), sizeof(uint64_t), 1);
+        CHECK(rc == 1);
+        CHECK(a == 1);
+
+        rc = mliIo_read(&buff, (const void *)(&b), sizeof(int32_t), 1);
+        CHECK(rc == 1);
+        CHECK(b == 42);
+
+        rc = mliIo_read(&buff, (const void *)(&b), sizeof(int32_t), 1);
+        CHECK(rc == 1);
+        CHECK(b == 13);
+
+        rc = mliIo_read(&buff, (const void *)(&c), sizeof(double), 1);
+        CHECK(rc == 1);
+        CHECK_MARGIN(c, 20 * 0.5, 1e-6);
+
+        mliIo_free(&buff);
+}

libs/mli_corsika/src/mli_corsika_version.h

@@ -4,6 +4,6 @@
 
 #define MLI_CORSIKA_VERSION_MAYOR 0
 #define MLI_CORSIKA_VERSION_MINOR 3
-#define MLI_CORSIKA_VERSION_PATCH 0
+#define MLI_CORSIKA_VERSION_PATCH 1
 
 #endif