shtr_line_view.c (11776B)
1 /* Copyright (C) 2022, 2025 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2025 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 18 19 #define _POSIX_C_SOURCE 200112L /* nextafter */ 20 21 #include "shtr.h" 22 #include "shtr_c.h" 23 #include "shtr_line_list_c.h" 24 25 #include <rsys/cstr.h> 26 #include <rsys/dynamic_array_size_t.h> 27 28 #include <math.h> 29 30 /* Version of the line view. One should increment it and perform a version 31 * management onto serialized data when the line view structure is updated. */ 32 static const int SHTR_LINE_VIEW_VERSION = 0; 33 34 struct shtr_line_view { 35 struct shtr_line_list* list; 36 struct darray_size_t line_ids; /* Indices of the selected lines */ 37 ref_T ref; 38 }; 39 40 struct molecule_selection { 41 /* Map the isotope local identifier to a boolean defining if the isotope is 42 * selected or not */ 43 char isotopes[SHTR_MAX_ISOTOPES_COUNT]; 44 double cutoff; /* Molecule cutoff */ 45 }; 46 47 /******************************************************************************* 48 * Helper functions 49 ******************************************************************************/ 50 static res_T 51 check_shtr_isotope_selection 52 (struct shtr* shtr, 53 const char* caller, 54 const struct shtr_isotope_selection* molecule) 55 { 56 size_t i; 57 ASSERT(caller && molecule); 58 59 if((size_t)molecule->id >= SHTR_MAX_MOLECULES_COUNT) { 60 ERROR(shtr, 61 "%s: molecule %d: invalid molecule identifier. " 62 "It must be less than %d.\n", 63 caller, molecule->id, SHTR_MAX_MOLECULES_COUNT); 64 return RES_BAD_ARG; 65 } 66 67 if(molecule->cutoff <= 0) { 68 ERROR(shtr, "%s: molecule %d: invalid cutoff %g.\n", 69 caller, molecule->id, molecule->cutoff); 70 return RES_BAD_ARG; 71 } 72 73 FOR_EACH(i, 0, molecule->nisotopes) { 74 if(molecule->isotope_ids_local[i] >= SHTR_MAX_ISOTOPES_COUNT) { 75 ERROR(shtr, 76 "%s: molecule %d: isotope %d: invalid isotope local identifier. " 77 "It must be less than %d.\n", 78 caller, 79 molecule->id, 80 molecule->isotope_ids_local[i], 81 SHTR_MAX_ISOTOPES_COUNT); 82 return RES_BAD_ARG; 83 } 84 } 85 return RES_OK; 86 } 87 88 89 static res_T 90 check_shtr_line_view_create_args 91 (struct shtr* shtr, 92 const char* caller, 93 const struct shtr_line_view_create_args* args) 94 { 95 size_t i; 96 ASSERT(caller); 97 98 if(!args) return RES_BAD_ARG; 99 100 if(args->wavenumber_range[0] > args->wavenumber_range[1]) { 101 ERROR(shtr, "%s: invalid line view spectral range [%g, %g].\n", 102 caller, args->wavenumber_range[0], args->wavenumber_range[1]); 103 return RES_BAD_ARG; 104 } 105 106 if(args->pressure < 0) { 107 ERROR(shtr, "%s: invalid pressure %g.\n", caller, args->pressure); 108 return RES_BAD_ARG; 109 } 110 111 FOR_EACH(i, 0, args->nmolecules) { 112 const res_T res = check_shtr_isotope_selection 113 (shtr, caller, &args->molecules[i]); 114 if(res != RES_OK) return res; 115 116 } 117 return RES_OK; 118 } 119 120 static res_T 121 create_line_view(struct shtr* shtr, struct shtr_line_view** out_view) 122 { 123 struct shtr_line_view* view = NULL; 124 res_T res = RES_OK; 125 ASSERT(shtr && out_view); 126 127 view = MEM_CALLOC(shtr->allocator, 1, sizeof(*view)); 128 if(!view) { 129 ERROR(shtr, "Could not allocate the line view.\n"); 130 res = RES_MEM_ERR; 131 goto error; 132 } 133 ref_init(&view->ref); 134 darray_size_t_init(shtr->allocator, &view->line_ids); 135 136 exit: 137 *out_view = view; 138 return res; 139 error: 140 if(view) { 141 SHTR(line_view_ref_put(view)); 142 view = NULL; 143 } 144 goto exit; 145 } 146 147 static res_T 148 select_lines 149 (struct shtr_line_view* view, 150 const char* caller, 151 const struct shtr_line_view_create_args* args) 152 { 153 const struct shtr_line* lines; 154 struct molecule_selection selection[SHTR_MAX_MOLECULES_COUNT]; 155 size_t imol; 156 size_t iiso; 157 size_t iline; 158 size_t nlines; 159 res_T res = RES_OK; 160 ASSERT(view && caller && args); 161 162 /* Nothing to do */ 163 if(args->nmolecules == 0) goto exit; 164 165 /* Setup the selection lookup table that map the isotope of a molecule to a 166 * boolean defining if it is selected */ 167 memset(selection, 0, sizeof(selection)); 168 FOR_EACH(imol, 0, args->nmolecules) { 169 const int32_t mol_id = args->molecules[imol].id; 170 ASSERT(mol_id < SHTR_MAX_MOLECULES_COUNT); 171 selection[mol_id].cutoff = args->molecules[imol].cutoff; 172 173 if(args->molecules[imol].nisotopes == 0) { /* Select all isotopes */ 174 FOR_EACH(iiso, 0, SHTR_MAX_ISOTOPES_COUNT) { 175 selection[mol_id].isotopes[iiso] = 1; 176 } 177 178 } else { 179 FOR_EACH(iiso, 0, args->molecules[imol].nisotopes) { 180 const int32_t iso_id = args->molecules[imol].isotope_ids_local[iiso]; 181 ASSERT(iso_id < SHTR_MAX_ISOTOPES_COUNT); 182 selection[mol_id].isotopes[iso_id] = 1; 183 } 184 } 185 } 186 187 lines = darray_line_cdata_get(&view->list->lines); 188 nlines = darray_line_size_get(&view->list->lines); 189 190 /* Iterate through list of lines to find the ones to select based on spectral 191 * range and isotope selection */ 192 FOR_EACH(iline, 0, nlines) { 193 const struct shtr_line* line = lines + iline; 194 double nu = 0; 195 196 /* The line is not selected */ 197 if(selection[line->molecule_id].isotopes[line->isotope_id_local] == 0) { 198 continue; 199 } 200 201 /* Compute the line center for the submitted pressure */ 202 nu = line->wavenumber + line->delta_air * args->pressure; 203 204 /* The line is out of the spectral range */ 205 if(nu + selection[line->molecule_id].cutoff < args->wavenumber_range[0] 206 || nu - selection[line->molecule_id].cutoff > args->wavenumber_range[1]) { 207 continue; 208 } 209 210 res = darray_size_t_push_back(&view->line_ids, &iline); 211 if(res != RES_OK) { 212 ERROR(view->list->shtr, 213 "%s: could not register the line into the view -- %s.\n", 214 caller, res_to_cstr(res)); 215 goto error; 216 } 217 } 218 219 exit: 220 return res; 221 error: 222 darray_size_t_clear(&view->line_ids); 223 goto exit; 224 } 225 226 static void 227 release_line_view(ref_T* ref) 228 { 229 struct shtr_line_view* view = CONTAINER_OF(ref, struct shtr_line_view, ref); 230 struct shtr_line_list* list = view->list; 231 ASSERT(ref); 232 darray_size_t_release(&view->line_ids); 233 MEM_RM(list->shtr->allocator, view); 234 SHTR(line_list_ref_put(list)); 235 } 236 237 /******************************************************************************* 238 * Exported functions 239 ******************************************************************************/ 240 res_T 241 shtr_line_view_create 242 (struct shtr_line_list* list, 243 const struct shtr_line_view_create_args* args, 244 struct shtr_line_view** out_view) 245 { 246 struct shtr_line_view* view = NULL; 247 res_T res = RES_OK; 248 249 if(!list || !out_view) { res = RES_BAD_ARG; goto error; } 250 251 res = check_shtr_line_view_create_args(list->shtr, FUNC_NAME, args); 252 if(res != RES_OK) goto error; 253 254 res = create_line_view(list->shtr, &view); 255 if(res != RES_OK) goto error; 256 SHTR(line_list_ref_get(list)); 257 view->list = list; 258 259 res = select_lines(view, FUNC_NAME, args); 260 if(res != RES_OK) goto error; 261 262 exit: 263 if(out_view) *out_view = view; 264 return res; 265 error: 266 if(view) { SHTR(line_view_ref_put(view)); view = NULL; } 267 goto exit; 268 } 269 270 res_T 271 shtr_line_view_create_from_stream 272 (struct shtr* shtr, 273 FILE* stream, 274 struct shtr_line_view** out_view) 275 { 276 struct shtr_line_view* view = NULL; 277 size_t nids; 278 int version; 279 res_T res = RES_OK; 280 281 if(!shtr || !stream || !out_view) { res = RES_BAD_ARG; goto error; } 282 283 res = create_line_view(shtr, &view); 284 if(res != RES_OK) goto error; 285 286 #define READ(Var, Nb) { \ 287 if(fread((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \ 288 if(feof(stream)) { \ 289 res = RES_BAD_ARG; \ 290 } else if(ferror(stream)) { \ 291 res = RES_IO_ERR; \ 292 } else { \ 293 res = RES_UNKNOWN_ERR; \ 294 } \ 295 ERROR(shtr, "%s: error reading isotope metadata -- %s.\n", \ 296 FUNC_NAME, res_to_cstr(res)); \ 297 goto error; \ 298 } \ 299 } (void)0 300 READ(&version, 1); 301 if(version != SHTR_LINE_VIEW_VERSION) { 302 ERROR(shtr, 303 "%s: unexpected line view version %d. " 304 "Expecting a line view in version %d.\n", 305 FUNC_NAME, version, SHTR_LINE_VIEW_VERSION); 306 res = RES_BAD_ARG; 307 goto error; 308 } 309 310 res = shtr_line_list_create_from_stream(shtr, stream, &view->list); 311 if(res != RES_OK) goto error; 312 313 READ(&nids, 1); 314 res = darray_size_t_resize(&view->line_ids, nids); 315 if(res != RES_OK) { 316 ERROR(shtr, 317 "%s: error allocating the line indices of the line view -- %s.\n", 318 FUNC_NAME, res_to_cstr(res)); 319 goto error; 320 } 321 322 READ(darray_size_t_data_get(&view->line_ids), nids); 323 #undef READ 324 325 exit: 326 if(out_view) *out_view = view; 327 return res; 328 error: 329 if(view) { 330 SHTR(line_view_ref_put(view)); 331 view = NULL; 332 } 333 goto exit; 334 } 335 336 res_T 337 shtr_line_view_ref_get(struct shtr_line_view* view) 338 { 339 if(!view) return RES_BAD_ARG; 340 ref_get(&view->ref); 341 return RES_OK; 342 } 343 344 res_T 345 shtr_line_view_ref_put(struct shtr_line_view* view) 346 { 347 if(!view) return RES_BAD_ARG; 348 ref_put(&view->ref, release_line_view); 349 return RES_OK; 350 } 351 352 res_T 353 shtr_line_view_get_size(const struct shtr_line_view* view, size_t* sz) 354 { 355 if(!view || !sz) return RES_BAD_ARG; 356 *sz = darray_size_t_size_get(&view->line_ids); 357 return RES_OK; 358 } 359 360 res_T 361 shtr_line_view_get_line 362 (const struct shtr_line_view* view, 363 const size_t iline, 364 const struct shtr_line** line) 365 { 366 size_t i; 367 if(!view || !line || iline >= darray_size_t_size_get(&view->line_ids)) { 368 return RES_BAD_ARG; 369 } 370 i = darray_size_t_cdata_get(&view->line_ids)[iline]; 371 *line = darray_line_cdata_get(&view->list->lines) + i; 372 return RES_OK; 373 } 374 375 res_T 376 shtr_line_view_write 377 (const struct shtr_line_view* view, 378 FILE* stream) 379 { 380 size_t nids = 0; 381 res_T res = RES_OK; 382 383 if(!view || !stream) { 384 res = RES_BAD_ARG; 385 goto error; 386 } 387 388 #define WRITE(Var, Nb) { \ 389 if(fwrite((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \ 390 ERROR(view->list->shtr, "%s: error writing line view.\n", FUNC_NAME); \ 391 res = RES_IO_ERR; \ 392 goto error; \ 393 } \ 394 } (void)0 395 WRITE(&SHTR_LINE_VIEW_VERSION, 1); 396 397 res = shtr_line_list_write(view->list, stream); 398 if(res != RES_OK) goto error; 399 400 nids = darray_size_t_size_get(&view->line_ids); 401 WRITE(&nids, 1); 402 WRITE(darray_size_t_cdata_get(&view->line_ids), nids); 403 #undef WRITE 404 405 exit: 406 return res; 407 error: 408 goto exit; 409 }