test_swf_HXd.h (4531B)
1 /* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "swf.h" 17 18 #include <rsys/rsys.h> 19 #include <rsys/math.h> 20 21 #include <stdlib.h> 22 23 #ifndef TEST_SWF_DIMENSION 24 #error "The TEST_SWF_DIMENSION macro must be defined" 25 #endif 26 27 /* 2D H function */ 28 #if TEST_SWF_DIMENSION == 2 29 #define SWF_HXD_TABULATE_ARGS_DEFAULT SWF_H2D_TABULATE_ARGS_DEFAULT 30 #define swf_HXd_eval swf_H2d_eval 31 #define swf_HXd_inverse swf_H2d_inverse 32 #define swf_HXd_tabulate swf_H2d_tabulate 33 34 /* 3D H function */ 35 #elif TEST_SWF_DIMENSION == 3 36 #define SWF_HXD_TABULATE_ARGS_DEFAULT SWF_H3D_TABULATE_ARGS_DEFAULT 37 #define swf_HXd_eval swf_H3d_eval 38 #define swf_HXd_inverse swf_H3d_inverse 39 #define swf_HXd_tabulate swf_H3d_tabulate 40 #else 41 #error "Invalid TEST_SWF_DIMENSION" 42 #endif 43 44 /* Macro defining a generic function name for the dimension*/ 45 #define Xd(Name) CONCAT(CONCAT(CONCAT(Name, _), TEST_SWF_DIMENSION), d) 46 47 static INLINE void 48 Xd(check_tabulation_creation)(void) 49 { 50 struct swf_H_tabulate_args args = SWF_HXD_TABULATE_ARGS_DEFAULT; 51 struct swf_tabulation* tab = NULL; 52 53 CHK(swf_HXd_tabulate(NULL, &tab) == RES_BAD_ARG); 54 CHK(swf_HXd_tabulate(&args, NULL) == RES_BAD_ARG); 55 CHK(swf_HXd_tabulate(&args, &tab) == RES_OK); 56 57 CHK(swf_tabulation_ref_get(NULL) == RES_BAD_ARG); 58 CHK(swf_tabulation_ref_get(tab) == RES_OK); 59 CHK(swf_tabulation_ref_put(NULL) == RES_BAD_ARG); 60 CHK(swf_tabulation_ref_put(tab) == RES_OK); 61 CHK(swf_tabulation_ref_put(tab) == RES_OK); 62 63 args.x_min = SWF_HXD_TABULATE_ARGS_DEFAULT.x_max; 64 args.x_max = SWF_HXD_TABULATE_ARGS_DEFAULT.x_min; 65 CHK(swf_HXd_tabulate(&args, &tab) == RES_BAD_ARG); 66 67 args.x_max = args.x_min; 68 CHK(swf_HXd_tabulate(&args, &tab) == RES_BAD_ARG); 69 70 args.x_min = SWF_HXD_TABULATE_ARGS_DEFAULT.x_min; 71 args.x_max = SWF_HXD_TABULATE_ARGS_DEFAULT.x_max; 72 args.step = 0; 73 CHK(swf_HXd_tabulate(&args, &tab) == RES_BAD_ARG); 74 75 args.step = 1.0e-2; 76 CHK(swf_HXd_tabulate(&args, &tab) == RES_OK); 77 CHK(swf_tabulation_ref_put(tab) == RES_OK); 78 } 79 80 static void 81 Xd(check_inversion)(void) 82 { 83 FILE* fp = NULL; 84 const double nsteps = 10; 85 86 struct swf_H_tabulate_args args = SWF_HXD_TABULATE_ARGS_DEFAULT; 87 struct swf_tabulation* tab = NULL; 88 double pHx = 0; 89 double x0 = 0; 90 size_t i = 0; 91 92 fp=fopen("H"STR(TEST_SWF_DIMENSION)"d.dat", "w"); 93 CHK(fp != NULL); 94 95 CHK(swf_HXd_tabulate(&args, &tab) == RES_OK); 96 97 for(x0 = args.x_min; x0 <= args.x_max; x0 += x0*args.step) { 98 const double x1 = MMIN(x0 + x0*args.step, args.x_max); 99 const double delta_x = (x1 - x0) / (double)nsteps; 100 101 for(i = 0; i < nsteps; ++i) { 102 const double x = x0 + (double)i*delta_x; 103 const double Hx = swf_HXd_eval(x); 104 const double xl = swf_tabulation_inverse(tab, SWF_LINEAR, Hx); 105 const double xq = swf_tabulation_inverse(tab, SWF_QUADRATIC, Hx); 106 const double xi = swf_HXd_inverse(Hx); 107 const double errl = fabs((x - xl) / x); 108 const double errq = fabs((x - xq) / x); 109 const double erri = fabs((x - xi) / x); 110 111 /* Do not check the value of x if the corresponding H is indistinguishable 112 * from an H calculated from a previous x. In other words, the H's are 113 * numerically equal and any of the corresponding x-values is a valid 114 * inversion result, numerically speaking. */ 115 if(Hx == pHx) continue; 116 117 fprintf(fp, "%e %e %e %e %e\n", x, Hx, errq, errl, erri); 118 if(1e-9 < Hx && Hx < (1.0 - 1e-9)) { 119 CHK(errl < 1.e-5); 120 CHK(errq < 1.e-7); 121 CHK(erri < 1.e-7); 122 } else { 123 CHK(errl < 1.e-3); 124 CHK(errq < 1.e-3); 125 CHK(erri < 1.e-3); 126 } 127 pHx = Hx; 128 } 129 } 130 131 CHK(fclose(fp) == 0); 132 CHK(swf_tabulation_ref_put(tab) == RES_OK); 133 } 134 135 /* Undefine macros used for genericity */ 136 #undef TEST_SWF_DIMENSION 137 #undef SWF_HXD_TABULATE_ARGS_DEFAULT 138 #undef swf_HXd_eval 139 #undef swf_HXd_inverse 140 #undef swf_HXd_tabulate 141 #undef Xd