From 8e6ffb2d84778ed0f70e05e98dd617bee1fd1f49 Mon Sep 17 00:00:00 2001 From: Sebastian Achim Mueller Date: Mon, 27 May 2024 09:50:16 +0200 Subject: [PATCH] make mliAvl more useable, turn ground_grid into a histogram. No longer store the bin assignment for each photon. --- libs/mli/src/mliAvlDict.c | 9 +- libs/mli/src/mliAvlDict.h | 1 + libs/mli/src/mliAvlTree.c | 1 + libs/mli/src/mli_version.h | 2 +- libs/mli_corsika/apps/ground_grid.main.c | 51 ++++++- .../mli_corsika/src/mli_corsika_Histogram2d.c | 144 ++++++++++++++++++ .../mli_corsika/src/mli_corsika_Histogram2d.h | 34 +++++ .../src/mli_corsika_Histogram2d.test.c | 58 +++++++ libs/mli_corsika/src/mli_corsika_version.h | 2 +- 9 files changed, 294 insertions(+), 8 deletions(-) create mode 100644 libs/mli_corsika/src/mli_corsika_Histogram2d.c create mode 100644 libs/mli_corsika/src/mli_corsika_Histogram2d.h create mode 100644 libs/mli_corsika/src/mli_corsika_Histogram2d.test.c diff --git a/libs/mli/src/mliAvlDict.c b/libs/mli/src/mliAvlDict.c index e99dda02..7862a174 100644 --- a/libs/mli/src/mliAvlDict.c +++ b/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; +} diff --git a/libs/mli/src/mliAvlDict.h b/libs/mli/src/mliAvlDict.h index e1f1332c..d1c06c74 100644 --- a/libs/mli/src/mliAvlDict.h +++ b/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 diff --git a/libs/mli/src/mliAvlTree.c b/libs/mli/src/mliAvlTree.c index 41e78380..f6d04bd2 100644 --- a/libs/mli/src/mliAvlTree.c +++ b/libs/mli/src/mliAvlTree.c @@ -1,4 +1,5 @@ #include "mliAvlTree.h" +#include /* Adopted from: * diff --git a/libs/mli/src/mli_version.h b/libs/mli/src/mli_version.h index 6bdc8de8..b0c88f3d 100644 --- a/libs/mli/src/mli_version.h +++ b/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); diff --git a/libs/mli_corsika/apps/ground_grid.main.c b/libs/mli_corsika/apps/ground_grid.main.c index e7f9dd74..dd0df3e8 100644 --- a/libs/mli_corsika/apps/ground_grid.main.c +++ b/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); diff --git a/libs/mli_corsika/src/mli_corsika_Histogram2d.c b/libs/mli_corsika/src/mli_corsika_Histogram2d.c new file mode 100644 index 00000000..c752a9d0 --- /dev/null +++ b/libs/mli_corsika/src/mli_corsika_Histogram2d.c @@ -0,0 +1,144 @@ +/* Copyright 2017 Sebastian A. Mueller */ + +#include "mli_corsika_Histogram2d.h" +#include +#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); +} diff --git a/libs/mli_corsika/src/mli_corsika_Histogram2d.h b/libs/mli_corsika/src/mli_corsika_Histogram2d.h new file mode 100644 index 00000000..5bfa34b2 --- /dev/null +++ b/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 +#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 diff --git a/libs/mli_corsika/src/mli_corsika_Histogram2d.test.c b/libs/mli_corsika/src/mli_corsika_Histogram2d.test.c new file mode 100644 index 00000000..20b9869c --- /dev/null +++ b/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); +} diff --git a/libs/mli_corsika/src/mli_corsika_version.h b/libs/mli_corsika/src/mli_corsika_version.h index d5fca9e5..3a6f6bc1 100644 --- a/libs/mli_corsika/src/mli_corsika_version.h +++ b/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