htcp

Properties of water suspended in clouds
git clone git://git.meso-star.fr/htcp.git
Log | Files | Refs | README | LICENSE

htcp.c (11868B)


      1 /* Copyright (C) 2018, 2020-2023, 2025, 2025 |Méso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2018 Centre National de la Recherche Scientifique
      3  * Copyright (C) 2018 Université Paul Sabatier
      4  *
      5  * This program is free software: you can redistribute it and/or modify
      6  * it under the terms of the GNU General Public License as published by
      7  * the Free Software Foundation, either version 3 of the License, or
      8  * (at your option) any later version.
      9  *
     10  * This program is distributed in the hope that it will be useful,
     11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     13  * GNU General Public License for more details.
     14  *
     15  * You should have received a copy of the GNU General Public License
     16  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     17 
     18 #define _POSIX_C_SOURCE 200809L /* mmap support */
     19 #define _DEFAULT_SOURCE 1 /* MAP_POPULATE support */
     20 #define _BSD_SOURCE 1 /* MAP_POPULATE for glibc < 2.19 */
     21 
     22 #include "htcp.h"
     23 
     24 #include <rsys/dynamic_array_double.h>
     25 #include <rsys/logger.h>
     26 #include <rsys/ref_count.h>
     27 #include <rsys/mem_allocator.h>
     28 
     29 #include <errno.h>
     30 #include <string.h>
     31 #include <sys/mman.h>
     32 #include <unistd.h>
     33 
     34 enum { X, Y, Z, TIME }; /* Helper constants */
     35 
     36 struct htcp {
     37   int64_t pagesize;
     38   int8_t irregular_z;
     39   int32_t definition[4];
     40   double lower[3];
     41   double upper[3];
     42   double vxsz[2]; /* Size of the voxels in X and Y */
     43   struct darray_double vxsz_z; /* Size of the voxels along the Z dimension */
     44   struct darray_double coord_z; /* Lower coordinate of the voxel along Z */
     45 
     46   double* RCT; /* Mapped memory */
     47   double* RVT; /* Mapped memory */
     48   double* PABST; /* Mapped memory */
     49   double* T; /* Mapped memory */
     50   size_t RCT_length; /* In bytes */
     51   size_t RVT_length; /* In bytes */
     52   size_t PABST_length; /* In bytes */
     53   size_t T_length; /* In bytes */
     54 
     55   size_t pagesize_os; /* Page size of the os */
     56   int verbose; /* Verbosity level */
     57   struct logger* logger;
     58   struct mem_allocator* allocator;
     59   ref_T ref;
     60 };
     61 
     62 /*******************************************************************************
     63  * Local functions
     64  ******************************************************************************/
     65 static void
     66 log_msg
     67   (const struct htcp* htcp,
     68    const enum log_type stream,
     69    const char* msg,
     70    va_list vargs)
     71 {
     72   ASSERT(htcp && msg);
     73   if(htcp->verbose) {
     74     res_T res; (void)res;
     75     res = logger_vprint(htcp->logger, stream, msg, vargs);
     76     ASSERT(res == RES_OK);
     77   }
     78 }
     79 
     80 static void
     81 log_err(const struct htcp* htcp, const char* msg, ...)
     82 {
     83   va_list vargs_list;
     84   ASSERT(htcp && msg);
     85   va_start(vargs_list, msg);
     86   log_msg(htcp, LOG_ERROR, msg, vargs_list);
     87   va_end(vargs_list);
     88 }
     89 
     90 static void
     91 reset_htcp(struct htcp* htcp)
     92 {
     93   ASSERT(htcp);
     94   htcp->pagesize = 0;
     95   htcp->irregular_z = 0;
     96   htcp->definition[0] = 0;
     97   htcp->definition[1] = 0;
     98   htcp->definition[2] = 0;
     99   htcp->lower[0] = -1;
    100   htcp->lower[1] = -1;
    101   htcp->lower[2] = -1;
    102   htcp->vxsz[0] = -1;
    103   htcp->vxsz[1] = -1;
    104   darray_double_clear(&htcp->vxsz_z);
    105   darray_double_clear(&htcp->coord_z);
    106   if(htcp->RCT) munmap(htcp->RCT, htcp->RCT_length);
    107   if(htcp->RVT) munmap(htcp->RVT, htcp->RVT_length);
    108   if(htcp->PABST) munmap(htcp->PABST, htcp->PABST_length);
    109   if(htcp->T) munmap(htcp->T, htcp->T_length);
    110   htcp->RCT = NULL;
    111   htcp->RVT = NULL;
    112   htcp->RCT_length = 0;
    113   htcp->RVT_length = 0;
    114 }
    115 
    116 static res_T
    117 load_stream(struct htcp* htcp, FILE* stream, const char* stream_name)
    118 {
    119   size_t nz = 0;
    120   size_t map_len = 0;
    121   size_t filesz;
    122   off_t offset = 0;
    123   res_T res = RES_OK;
    124   ASSERT(htcp && stream && stream_name);
    125 
    126   reset_htcp(htcp);
    127 
    128   #define READ(Var, N, Name) {                                                 \
    129     if(fread((Var), sizeof(*(Var)), (N), stream) != (N)) {                     \
    130       log_err(htcp, "%s: could not read the %s\n", stream_name, Name);         \
    131       res = RES_IO_ERR;                                                        \
    132       goto error;                                                              \
    133     }                                                                          \
    134   } (void)0
    135   READ(&htcp->pagesize, 1, "page size");
    136   if(!IS_POW2(htcp->pagesize)) {
    137     log_err(htcp, "%s: invalid page size `%li'. It must be a power of two.\n",
    138       stream_name, (long)htcp->pagesize);
    139     res = RES_BAD_ARG;
    140     goto error;
    141   }
    142 
    143   if(!IS_ALIGNED(htcp->pagesize, htcp->pagesize_os)) {
    144     log_err(htcp,
    145 "%s: invalid page size `%li'. The page size attribute must be aligned on the "
    146 "page size of the operating system , i.e. %lu.\n",
    147       stream_name, (long)htcp->pagesize,  (unsigned long)htcp->pagesize_os);
    148     res = RES_BAD_ARG;
    149     goto error;
    150   }
    151 
    152   READ(&htcp->irregular_z, 1, "'irregular Z' flag");
    153   READ(htcp->definition, 4, "spatial and time definitions");
    154   if(htcp->definition[0] <= 0
    155   || htcp->definition[1] <= 0
    156   || htcp->definition[2] <= 0
    157   || htcp->definition[3] <= 0) {
    158     log_err(htcp,
    159 "%s: the spatial/time definition cannot be negative or null -- spatial "
    160 "definition: %i %i %i; time definition: %i\n",
    161       stream_name, SPLIT3(htcp->definition), htcp->definition[TIME]);
    162     res = RES_BAD_ARG;
    163     goto error;
    164   }
    165 
    166   READ(htcp->lower, 3, "lower position");
    167   READ(htcp->vxsz, 2, "XY voxel size ");
    168 
    169   nz = htcp->irregular_z ? (size_t)htcp->definition[Z] : 1;
    170   res = darray_double_resize(&htcp->vxsz_z, nz);
    171   if(res != RES_OK) {
    172     log_err(htcp,
    173       "%s: could not allocate memory to store the size of the voxels in Z.\n",
    174       stream_name);
    175     goto error;
    176   }
    177   READ(darray_double_data_get(&htcp->vxsz_z), nz, "Z voxel size(s)");
    178   #undef READ
    179 
    180   htcp->upper[0] = htcp->lower[0] + htcp->vxsz[0] * htcp->definition[0];
    181   htcp->upper[1] = htcp->lower[1] + htcp->vxsz[1] * htcp->definition[1];
    182   if(!htcp->irregular_z) {
    183     htcp->upper[2] = htcp->lower[2]
    184       + darray_double_cdata_get(&htcp->vxsz_z)[0] * htcp->definition[2];
    185   } else {
    186     /* Compute the Z lower bound in Z of each Z slice */
    187     const double* size = NULL;
    188     double* coord = NULL;
    189     size_t i;
    190     res = darray_double_resize(&htcp->coord_z, nz);
    191     if(res != RES_OK) {
    192       log_err(htcp, "%s: could not allocate memory to store the lower "
    193         "coordinate of the voxels in Z.\n", FUNC_NAME);
    194       goto error;
    195     }
    196     size  = darray_double_cdata_get(&htcp->vxsz_z);
    197     coord = darray_double_data_get(&htcp->coord_z);
    198     FOR_EACH(i, 0, nz) coord[i] = i ? coord[i-1] + size[i-1] : htcp->lower[2];
    199     htcp->upper[2] = coord[nz-1] + size[nz-1];
    200   }
    201 
    202   map_len =
    203     (size_t)htcp->definition[X]
    204   * (size_t)htcp->definition[Y]
    205   * (size_t)htcp->definition[Z]
    206   * (size_t)htcp->definition[TIME]
    207   * sizeof(double);
    208   /* Align sizeof mapped data onto page size */
    209   map_len = ALIGN_SIZE(map_len, (size_t)htcp->pagesize);
    210 
    211   /* Ensure that offset is align on page size */
    212   offset = ALIGN_SIZE(ftell(stream), htcp->pagesize);
    213 
    214   fseek(stream, 0, SEEK_END);
    215   filesz = (size_t)ftell(stream);
    216 
    217   #define MMAP(Var) {                                                          \
    218     if((size_t)offset + map_len > filesz) {                                    \
    219       log_err(htcp, "%s: coult not load the "STR(Var)" data.\n", stream_name); \
    220       res = RES_IO_ERR;                                                        \
    221       goto error;                                                              \
    222     }                                                                          \
    223     htcp->Var = mmap(NULL, map_len, PROT_READ, MAP_PRIVATE|MAP_POPULATE,       \
    224       fileno(stream), offset);                                                 \
    225     if(htcp->Var == MAP_FAILED) {                                              \
    226       log_err(htcp, "%s: could not map the "STR(Var)" data -- %s.\n",          \
    227         stream_name, strerror(errno));                                         \
    228       res = RES_IO_ERR;                                                        \
    229       goto error;                                                              \
    230     }                                                                          \
    231     htcp->Var##_length = map_len;                                              \
    232   } (void)0
    233   MMAP(RVT); offset += (off_t)map_len;
    234   MMAP(RCT); offset += (off_t)map_len;
    235   MMAP(PABST); offset += (off_t)map_len;
    236   MMAP(T); offset += (off_t)map_len;
    237   #undef MMAP
    238 
    239   CHK(fseek(stream, offset, SEEK_CUR) != -1);
    240 
    241 exit:
    242   return res;
    243 error:
    244   goto exit;
    245 }
    246 
    247 static void
    248 release_htcp(ref_T* ref)
    249 {
    250   struct htcp* htcp;
    251   ASSERT(ref);
    252   htcp = CONTAINER_OF(ref, struct htcp, ref);
    253   reset_htcp(htcp);
    254   darray_double_release(&htcp->vxsz_z);
    255   darray_double_release(&htcp->coord_z);
    256   MEM_RM(htcp->allocator, htcp);
    257 }
    258 
    259 /*******************************************************************************
    260  * Exported functions
    261  ******************************************************************************/
    262 res_T
    263 htcp_create
    264   (struct logger* in_logger,
    265    struct mem_allocator* mem_allocator,
    266    const int verbose,
    267    struct htcp** out_htcp)
    268 {
    269   struct htcp* htcp = NULL;
    270   struct mem_allocator* allocator = NULL;
    271   struct logger* logger = NULL;
    272   res_T res = RES_OK;
    273 
    274   if(!out_htcp) {
    275     res = RES_BAD_ARG;
    276     goto error;
    277   }
    278 
    279   allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
    280   logger = in_logger ? in_logger : LOGGER_DEFAULT;
    281 
    282   htcp = MEM_CALLOC(allocator, 1, sizeof(struct htcp));
    283   if(!htcp) {
    284     if(verbose) {
    285       logger_print(logger, LOG_ERROR,
    286         "%s: could not allocate the HTCP handler.\n", FUNC_NAME);
    287     }
    288     res = RES_MEM_ERR;
    289     goto error;
    290   }
    291   ref_init(&htcp->ref);
    292   htcp->allocator = allocator;
    293   htcp->logger = logger;
    294   htcp->verbose = verbose;
    295   htcp->pagesize_os = (size_t)sysconf(_SC_PAGESIZE);
    296   darray_double_init(htcp->allocator, &htcp->vxsz_z);
    297   darray_double_init(htcp->allocator, &htcp->coord_z);
    298 
    299 exit:
    300   if(out_htcp) *out_htcp = htcp;
    301   return res;
    302 error:
    303   if(htcp) {
    304     HTCP(ref_put(htcp));
    305     htcp = NULL;
    306   }
    307   goto exit;
    308 }
    309 
    310 res_T
    311 htcp_ref_get(struct htcp* htcp)
    312 {
    313   if(!htcp) return RES_BAD_ARG;
    314   ref_get(&htcp->ref);
    315   return RES_OK;
    316 }
    317 
    318 res_T
    319 htcp_ref_put(struct htcp* htcp)
    320 {
    321   if(!htcp) return RES_BAD_ARG;
    322   ref_put(&htcp->ref, release_htcp);
    323   return RES_OK;
    324 }
    325 
    326 res_T
    327 htcp_load(struct htcp* htcp, const char* path)
    328 {
    329   FILE* file = NULL;
    330   res_T res = RES_OK;
    331 
    332   if(!htcp || !path) {
    333     res = RES_BAD_ARG;
    334     goto error;
    335   }
    336 
    337   file = fopen(path, "r");
    338   if(!file) {
    339     log_err(htcp, "%s: error opening file `%s'.\n", FUNC_NAME, path);
    340     res = RES_IO_ERR;
    341     goto error;
    342   }
    343 
    344   res = load_stream(htcp, file, path);
    345   if(res != RES_OK) goto error;
    346 
    347 exit:
    348   if(file) fclose(file);
    349   return res;
    350 error:
    351   goto exit;
    352 }
    353 
    354 res_T
    355 htcp_load_stream(struct htcp* htcp, FILE* stream)
    356 {
    357   if(!htcp || !stream) return RES_BAD_ARG;
    358   return load_stream(htcp, stream, "<stream>");
    359 }
    360 
    361 res_T
    362 htcp_get_desc(const struct htcp* htcp, struct htcp_desc* desc)
    363 {
    364   if(!htcp || !desc) return RES_BAD_ARG;
    365   if(!htcp->RVT) return RES_BAD_ARG;
    366   desc->pagesize = (size_t)htcp->pagesize;
    367   desc->irregular_z = htcp->irregular_z != 0;
    368   desc->spatial_definition[0] = (size_t)htcp->definition[X];
    369   desc->spatial_definition[1] = (size_t)htcp->definition[Y];
    370   desc->spatial_definition[2] = (size_t)htcp->definition[Z];
    371   desc->time_definition = (size_t)htcp->definition[TIME];
    372   desc->lower[0] = htcp->lower[0];
    373   desc->lower[1] = htcp->lower[1];
    374   desc->lower[2] = htcp->lower[2];
    375   desc->upper[0] = htcp->upper[0];
    376   desc->upper[1] = htcp->upper[1];
    377   desc->upper[2] = htcp->upper[2];
    378   desc->vxsz_x = htcp->vxsz[0];
    379   desc->vxsz_y = htcp->vxsz[1];
    380   desc->vxsz_z = darray_double_cdata_get(&htcp->vxsz_z);
    381   desc->coord_z = htcp->irregular_z ? darray_double_cdata_get(&htcp->coord_z) : NULL;
    382   desc->RCT = htcp->RCT;
    383   desc->RVT = htcp->RVT;
    384   desc->PABST = htcp->PABST;
    385   desc->T = htcp->T;
    386   return RES_OK;
    387 }
    388