test_shtr_isotope_metadata.c (15461B)
1 /* Copyright (C) 2022, 2025, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2025, 2026 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 #include "shtr.h" 20 21 #include <rsys/mem_allocator.h> 22 #include <rsys/math.h> 23 24 #include <string.h> 25 26 static const struct shtr_isotope H2O_isotopes[] = { 27 {9.97317E-01, 1.7458E+02, 18.010565, 0, 1, 161}, 28 {1.99983E-03, 1.7605E+02, 20.014811, 0, 1, 181}, 29 {3.71884E-04, 1.0521E+03, 19.014780, 0, 6, 171}, 30 {3.10693E-04, 8.6474E+02, 19.016740, 0, 6, 162}, 31 {6.23003E-07, 8.7557E+02, 21.020985, 0, 6, 182}, 32 {1.15853E-07, 5.2268E+03, 20.020956, 0, 36, 172}, 33 {2.41974E-08, 1.0278E+03, 20.022915, 0, 1, 262} 34 }; 35 36 37 static const struct shtr_molecule H2O = { 38 "H2O", sizeof(H2O_isotopes)/sizeof(struct shtr_isotope), H2O_isotopes, 1, 39 }; 40 41 static const struct shtr_isotope CO2_isotopes[] = { 42 {9.84204E-01, 2.8609E+02, 43.989830, 1, 1, 626}, 43 {1.10574E-02, 5.7664E+02, 44.993185, 1, 2, 636}, 44 {3.94707E-03, 6.0781E+02, 45.994076, 1, 1, 628}, 45 {7.33989E-04, 3.5426E+03, 44.994045, 1, 6, 627}, 46 {4.43446E-05, 1.2255E+03, 46.997431, 1, 2, 638}, 47 {8.24623E-06, 7.1413E+03, 45.997400, 1, 12, 637}, 48 {3.95734E-06, 3.2342E+02, 47.998320, 1, 1, 828}, 49 {1.47180E-06, 3.7666E+03, 46.998291, 1, 6, 827}, 50 {1.36847E-07, 1.0972E+04, 45.998262, 1, 1, 727}, 51 {4.44600E-08, 6.5224E+02, 49.001675, 1, 2, 838}, 52 {1.65354E-08, 7.5950E+03, 48.001646, 1, 12, 837}, 53 {1.53745E-09, 2.2120E+04, 47.001618, 1, 2, 737} 54 }; 55 56 static const struct shtr_molecule CO2 = { 57 "CO2", sizeof(CO2_isotopes)/sizeof(struct shtr_isotope), CO2_isotopes, 2, 58 }; 59 60 static void 61 isotope_print(FILE* fp, const struct shtr_isotope* isotope) 62 { 63 CHK(fp && isotope); 64 fprintf(fp, " %d %.5E %.4E %d %.6f\n", 65 isotope->id, 66 isotope->abundance, 67 isotope->Q296K, 68 isotope->gj, 69 isotope->molar_mass); 70 } 71 72 static void 73 molecule_print(FILE* fp, const struct shtr_molecule* molecule) 74 { 75 size_t i; 76 CHK(fp && molecule); 77 78 fprintf(fp, " %s (%d)\n", molecule->name, molecule->id); 79 FOR_EACH(i, 0, molecule->nisotopes) { 80 isotope_print(fp, molecule->isotopes+i); 81 } 82 } 83 84 static int 85 isotope_eq(const struct shtr_isotope* i0, const struct shtr_isotope* i1) 86 { 87 CHK(i0 && i1); 88 return i0->abundance == i1->abundance 89 && i0->Q296K == i1->Q296K 90 && i0->molar_mass == i1->molar_mass 91 && i0->molecule_id_local == i1->molecule_id_local 92 && i0->gj == i1->gj 93 && i0->id == i1->id; 94 } 95 96 static int 97 molecule_eq(const struct shtr_molecule* m0, const struct shtr_molecule* m1) 98 { 99 size_t i; 100 101 CHK(m0 && m1); 102 if(strcmp(m0->name, m1->name) 103 || m0->id != m1->id 104 || m0->nisotopes != m1->nisotopes) 105 return 0; 106 107 FOR_EACH(i, 0, m0->nisotopes) { 108 if(!isotope_eq(m0->isotopes+i, m1->isotopes+i)) 109 return 0; 110 } 111 return 1; 112 } 113 114 static void 115 check_isotope 116 (struct shtr_isotope_metadata* mdata, 117 const struct shtr_molecule* molecule, 118 const struct shtr_isotope* isotope) 119 { 120 struct shtr_molecule molecule2 = SHTR_MOLECULE_NULL; 121 122 /* Check NaN */ 123 CHK(isotope->abundance == isotope->abundance); 124 CHK(isotope->Q296K == isotope->Q296K); 125 CHK(isotope->molar_mass == isotope->molar_mass); 126 127 CHK(isotope->abundance > 0 && isotope->abundance <= 1); 128 CHK(isotope->Q296K > 0); 129 CHK(isotope->molar_mass > 0); 130 CHK(isotope->id >= 0); 131 132 CHK(shtr_isotope_metadata_get_molecule 133 (mdata, isotope->molecule_id_local, &molecule2) == RES_OK); 134 CHK(molecule_eq(molecule, &molecule2)); 135 } 136 137 static void 138 check_molecule 139 (struct shtr_isotope_metadata* mdata, 140 struct shtr_molecule* molecule) 141 { 142 size_t i = 0; 143 144 CHK(mdata && molecule); 145 CHK(molecule->name); 146 CHK(molecule->id >= 0); 147 148 FOR_EACH(i, 0, molecule->nisotopes) { 149 check_isotope(mdata, molecule, molecule->isotopes+i); 150 } 151 } 152 153 static void 154 test_load(struct shtr* shtr) 155 { 156 struct shtr_molecule molecule = SHTR_MOLECULE_NULL; 157 const char* filename = "test_isotope_metadata.txt"; 158 struct shtr_isotope_metadata* mdata = NULL; 159 FILE* fp = NULL; 160 size_t nmolecules = 0; 161 size_t nisotopes = 0; 162 163 CHK(fp = fopen(filename, "w+")); 164 165 fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); 166 molecule_print(fp, &H2O); 167 molecule_print(fp, &CO2); 168 rewind(fp); 169 170 CHK(shtr_isotope_metadata_load_stream(NULL, fp, NULL, &mdata) == RES_BAD_ARG); 171 CHK(shtr_isotope_metadata_load_stream(shtr, NULL, NULL, &mdata) == RES_BAD_ARG); 172 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, NULL) == RES_BAD_ARG); 173 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_OK); 174 175 CHK(shtr_isotope_metadata_get_molecules_count(NULL, &nmolecules) == RES_BAD_ARG); 176 CHK(shtr_isotope_metadata_get_molecules_count(mdata, NULL) == RES_BAD_ARG); 177 CHK(shtr_isotope_metadata_get_molecules_count(mdata, &nmolecules) == RES_OK); 178 CHK(nmolecules == 2); 179 180 CHK(shtr_isotope_metadata_get_isotopes_count(NULL, &nisotopes) == RES_BAD_ARG); 181 CHK(shtr_isotope_metadata_get_isotopes_count(mdata, NULL) == RES_BAD_ARG); 182 CHK(shtr_isotope_metadata_get_isotopes_count(mdata, &nisotopes) == RES_OK); 183 CHK(nisotopes == H2O.nisotopes + CO2.nisotopes); 184 185 CHK(shtr_isotope_metadata_get_molecule(NULL, 0, &molecule) == RES_BAD_ARG); 186 CHK(shtr_isotope_metadata_get_molecule(mdata, 2, &molecule) == RES_BAD_ARG); 187 188 CHK(shtr_isotope_metadata_get_molecule(mdata, 0, &molecule) == RES_OK); 189 CHK(molecule_eq(&molecule, &H2O)); 190 check_molecule(mdata, &molecule); 191 192 CHK(shtr_isotope_metadata_get_molecule(mdata, 1, &molecule) == RES_OK); 193 CHK(molecule_eq(&molecule, &CO2)); 194 check_molecule(mdata, &molecule); 195 196 CHK(shtr_isotope_metadata_find_molecule(NULL, 1, &molecule) == RES_BAD_ARG); 197 CHK(shtr_isotope_metadata_find_molecule(mdata, 1, NULL) == RES_BAD_ARG); 198 CHK(shtr_isotope_metadata_find_molecule(mdata, 1, &molecule) == RES_OK); 199 CHK(!SHTR_MOLECULE_IS_NULL(&molecule)); 200 CHK(molecule_eq(&molecule, &H2O)); 201 202 CHK(shtr_isotope_metadata_find_molecule(mdata, 2, &molecule) == RES_OK); 203 CHK(!SHTR_MOLECULE_IS_NULL(&molecule)); 204 CHK(molecule_eq(&molecule, &CO2)); 205 206 CHK(shtr_isotope_metadata_find_molecule(mdata, 0, &molecule) == RES_OK); 207 CHK(SHTR_MOLECULE_IS_NULL(&molecule)); 208 209 CHK(shtr_isotope_metadata_ref_get(NULL) == RES_BAD_ARG); 210 CHK(shtr_isotope_metadata_ref_get(mdata) == RES_OK); 211 CHK(shtr_isotope_metadata_ref_put(NULL) == RES_BAD_ARG); 212 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 213 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 214 215 CHK(fclose(fp) == 0); 216 217 CHK(shtr_isotope_metadata_load(NULL, filename, &mdata) == RES_BAD_ARG); 218 CHK(shtr_isotope_metadata_load(shtr, NULL, &mdata) == RES_BAD_ARG); 219 CHK(shtr_isotope_metadata_load(shtr, filename, NULL) == RES_BAD_ARG); 220 CHK(shtr_isotope_metadata_load(shtr, "no_file", &mdata) == RES_IO_ERR); 221 222 CHK(shtr_isotope_metadata_load(shtr, filename, &mdata) == RES_OK); 223 CHK(shtr_isotope_metadata_get_molecules_count(mdata, &nmolecules) == RES_OK); 224 CHK(nmolecules == 2); 225 CHK(shtr_isotope_metadata_get_molecule(mdata, 0, &molecule) == RES_OK); 226 CHK(molecule_eq(&molecule, &H2O)); 227 CHK(shtr_isotope_metadata_get_molecule(mdata, 1, &molecule) == RES_OK); 228 CHK(molecule_eq(&molecule, &CO2)); 229 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 230 231 /* Empty file */ 232 CHK(fp = tmpfile()); 233 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_OK); 234 CHK(shtr_isotope_metadata_get_molecules_count(mdata, &nmolecules) == RES_OK); 235 CHK(nmolecules == 0); 236 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 237 fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); 238 rewind(fp); 239 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_OK); 240 CHK(shtr_isotope_metadata_get_molecules_count(mdata, &nmolecules) == RES_OK); 241 CHK(nmolecules == 0); 242 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 243 CHK(fclose(fp) == 0); 244 245 /* Molecule without isotope */ 246 CHK(fp = tmpfile()); 247 fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); 248 fprintf(fp, "%s (%d)\n", H2O.name, H2O.id); 249 fprintf(fp, "%s (%d)\n", CO2.name, CO2.id); 250 rewind(fp); 251 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_OK); 252 CHK(shtr_isotope_metadata_get_molecules_count(mdata, &nmolecules) == RES_OK); 253 CHK(nmolecules == 2); 254 255 CHK(shtr_isotope_metadata_get_molecule(mdata, 0, &molecule) == RES_OK); 256 CHK(!strcmp(molecule.name, H2O.name)); 257 CHK(molecule.id == H2O.id); 258 CHK(molecule.nisotopes == 0); 259 CHK(molecule.isotopes == NULL); 260 261 CHK(shtr_isotope_metadata_get_molecule(mdata, 1, &molecule) == RES_OK); 262 CHK(!strcmp(molecule.name, CO2.name)); 263 CHK(molecule.id == CO2.id); 264 CHK(molecule.nisotopes == 0); 265 CHK(molecule.isotopes == NULL); 266 267 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 268 CHK(fclose(fp) == 0); 269 } 270 271 static void 272 check_equality 273 (const struct shtr_isotope_metadata* mdata1, 274 const struct shtr_isotope_metadata* mdata2) 275 { 276 size_t n1, n2; 277 size_t imol, nmol; 278 CHK(mdata1 && mdata2); 279 280 CHK(shtr_isotope_metadata_get_molecules_count(mdata1, &n1) == RES_OK); 281 CHK(shtr_isotope_metadata_get_molecules_count(mdata2, &n2) == RES_OK); 282 CHK(n1 == n2); 283 nmol = n1; 284 285 CHK(shtr_isotope_metadata_get_isotopes_count(mdata1, &n1) == RES_OK); 286 CHK(shtr_isotope_metadata_get_isotopes_count(mdata2, &n2) == RES_OK); 287 CHK(n1 == n2); 288 289 FOR_EACH(imol, 0, nmol) { 290 struct shtr_molecule mol1 = SHTR_MOLECULE_NULL; 291 struct shtr_molecule mol2 = SHTR_MOLECULE_NULL; 292 size_t iiso, niso; 293 294 CHK(shtr_isotope_metadata_get_molecule(mdata1, imol, &mol1) == RES_OK); 295 CHK(shtr_isotope_metadata_get_molecule(mdata2, imol, &mol2) == RES_OK); 296 CHK(!strcmp(mol1.name, mol2.name)); 297 CHK(mol1.id == mol2.id); 298 CHK(mol1.nisotopes == mol2.nisotopes); 299 niso = mol1.nisotopes; 300 301 FOR_EACH(iiso, 0, niso) { 302 CHK(mol1.isotopes[iiso].abundance == mol2.isotopes[iiso].abundance); 303 CHK(mol1.isotopes[iiso].Q296K == mol2.isotopes[iiso].Q296K); 304 CHK(mol1.isotopes[iiso].molar_mass == mol2.isotopes[iiso].molar_mass); 305 CHK(mol1.isotopes[iiso].molecule_id_local == mol2.isotopes[iiso].molecule_id_local); 306 CHK(mol1.isotopes[iiso].gj == mol2.isotopes[iiso].gj); 307 CHK(mol1.isotopes[iiso].id == mol2.isotopes[iiso].id); 308 } 309 } 310 } 311 312 static void 313 test_serialization(struct shtr* shtr) 314 { 315 struct shtr_isotope_metadata* mdata1 = NULL; 316 struct shtr_isotope_metadata* mdata2 = NULL; 317 FILE* fp = NULL; 318 319 CHK(fp = tmpfile()); 320 321 fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); 322 molecule_print(fp, &H2O); 323 molecule_print(fp, &CO2); 324 rewind(fp); 325 326 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata1) == RES_OK); 327 fclose(fp); 328 329 CHK(fp = tmpfile()); 330 CHK(shtr_isotope_metadata_write(NULL, fp) == RES_BAD_ARG); 331 CHK(shtr_isotope_metadata_write(mdata1, NULL) == RES_BAD_ARG); 332 CHK(shtr_isotope_metadata_write(mdata1, fp) == RES_OK); 333 rewind(fp); 334 335 CHK(shtr_isotope_metadata_create_from_stream(NULL, fp, &mdata2) == RES_BAD_ARG); 336 CHK(shtr_isotope_metadata_create_from_stream(shtr, NULL, &mdata2) == RES_BAD_ARG); 337 CHK(shtr_isotope_metadata_create_from_stream(shtr, fp, NULL) == RES_BAD_ARG); 338 CHK(shtr_isotope_metadata_create_from_stream(shtr, fp, &mdata2) == RES_OK); 339 fclose(fp); 340 341 check_equality(mdata1, mdata2); 342 343 CHK(shtr_isotope_metadata_ref_put(mdata1) == RES_OK); 344 CHK(shtr_isotope_metadata_ref_put(mdata2) == RES_OK); 345 } 346 347 static void 348 test_load_failures(struct shtr* shtr) 349 { 350 struct shtr_isotope isotope = SHTR_ISOTOPE_NULL; 351 struct shtr_isotope_metadata* mdata = NULL; 352 FILE* fp = NULL; 353 354 CHK(shtr); 355 356 /* 1st line is missing */ 357 CHK(fp = tmpfile()); 358 molecule_print(fp, &H2O); 359 rewind(fp); 360 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 361 CHK(fclose(fp) == 0); 362 363 /* Invalid molecule id */ 364 CHK(fp = tmpfile()); 365 fprintf(fp, "Comment line\n"); 366 fprintf(fp, "H2O 1\n"); 367 isotope_print(fp, &H2O_isotopes[0]); 368 rewind(fp); 369 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 370 rewind(fp); 371 fprintf(fp, "Comment line\n"); 372 fprintf(fp, "H2O (-1)\n"); 373 isotope_print(fp, &H2O_isotopes[0]); 374 rewind(fp); 375 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 376 CHK(fclose(fp) == 0); 377 378 /* Invalid isotope id */ 379 CHK(fp = tmpfile()); 380 fprintf(fp, "Comment line\n"); 381 fprintf(fp, "H2O (1)\n"); 382 isotope = H2O_isotopes[0]; 383 isotope.id = -1; 384 isotope_print(fp, &isotope); 385 rewind(fp); 386 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 387 CHK(fclose(fp) == 0); 388 389 /* Invalid abundance */ 390 CHK(fp = tmpfile()); 391 fprintf(fp, "Comment line\n"); 392 fprintf(fp, "H2O (1)\n"); 393 isotope = H2O_isotopes[0]; 394 isotope.abundance = 0; 395 isotope_print(fp, &isotope); 396 rewind(fp); 397 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 398 CHK(fclose(fp) == 0); 399 CHK(fp = tmpfile()); 400 fprintf(fp, "Comment line\n"); 401 fprintf(fp, "H2O (1)\n"); 402 isotope.abundance = 1.00001; 403 isotope_print(fp, &isotope); 404 rewind(fp); 405 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 406 CHK(fclose(fp) == 0); 407 408 /* Invalid Q(296K) */ 409 CHK(fp = tmpfile()); 410 fprintf(fp, "Comment line\n"); 411 fprintf(fp, "H2O (1)\n"); 412 isotope = H2O_isotopes[0]; 413 isotope.Q296K = 0; 414 isotope_print(fp, &isotope); 415 rewind(fp); 416 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 417 CHK(fclose(fp) == 0); 418 419 /* Invalid molar mass */ 420 CHK(fp = tmpfile()); 421 fprintf(fp, "Comment line\n"); 422 fprintf(fp, "H2O (1)\n"); 423 isotope = H2O_isotopes[0]; 424 isotope.molar_mass = 0; 425 isotope_print(fp, &isotope); 426 rewind(fp); 427 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &mdata) == RES_BAD_ARG); 428 CHK(fclose(fp) == 0); 429 } 430 431 static void 432 test_load_file(struct shtr* shtr, const char* path) 433 { 434 struct shtr_isotope_metadata* mdata = NULL; 435 size_t i, n; 436 CHK(path); 437 printf("Loading `%s'.\n", path); 438 CHK(shtr_isotope_metadata_load(shtr, path, &mdata) == RES_OK); 439 CHK(shtr_isotope_metadata_get_molecules_count(mdata, &n) == RES_OK); 440 printf(" #molecules: %lu\n", n); 441 FOR_EACH(i, 0, n) { 442 struct shtr_molecule molecule = SHTR_MOLECULE_NULL; 443 CHK(shtr_isotope_metadata_get_molecule(mdata, i, &molecule) == RES_OK); 444 printf(" Checking %s\n", molecule.name); 445 check_molecule(mdata, &molecule); 446 } 447 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 448 } 449 450 int 451 main(int argc, char** argv) 452 { 453 struct shtr_create_args args = SHTR_CREATE_ARGS_DEFAULT; 454 struct shtr* shtr = NULL; 455 int i; 456 (void)argc, (void)argv; 457 458 args.verbose = 1; 459 CHK(shtr_create(&args, &shtr) == RES_OK); 460 461 test_load(shtr); 462 test_load_failures(shtr); 463 test_serialization(shtr); 464 FOR_EACH(i, 1, argc) { 465 test_load_file(shtr, argv[i]); 466 } 467 468 CHK(shtr_ref_put(shtr) == RES_OK); 469 CHK(mem_allocated_size() == 0); 470 return 0; 471 }