test.c (9338B)
1 #include "../squiggle.h" 2 #include <math.h> 3 #include <stdint.h> 4 #include <stdio.h> 5 #include <stdlib.h> 6 7 #define TOLERANCE 5.0 / 1000.0 8 #define MAX_NAME_LENGTH 500 9 10 // Structs 11 12 struct array_expectations { 13 double* array; 14 int n; 15 char* name; 16 double expected_mean; 17 double expected_std; 18 double tolerance; 19 }; 20 21 void test_array_expectations(struct array_expectations e) 22 { 23 double mean = array_mean(e.array, e.n); 24 double delta_mean = mean - e.expected_mean; 25 26 double std = array_std(e.array, e.n); 27 double delta_std = std - e.expected_std; 28 29 if ((fabs(delta_mean) / fabs(mean) > e.tolerance) && (fabs(delta_mean) > e.tolerance)) { 30 printf("[-] Mean test for %s NOT passed.\n", e.name); 31 printf("Mean of %s: %f, vs expected mean: %f\n", e.name, mean, e.expected_mean); 32 printf("delta: %f, relative delta: %f\n", delta_mean, delta_mean / fabs(mean)); 33 } else { 34 printf("[x] Mean test for %s PASSED\n", e.name); 35 } 36 37 if ((fabs(delta_std) / fabs(std) > e.tolerance) && (fabs(delta_std) > e.tolerance)) { 38 printf("[-] Std test for %s NOT passed.\n", e.name); 39 printf("Std of %s: %f, vs expected std: %f\n", e.name, std, e.expected_std); 40 printf("delta: %f, relative delta: %f\n", delta_std, delta_std / fabs(std)); 41 } else { 42 printf("[x] Std test for %s PASSED\n", e.name); 43 } 44 45 printf("\n"); 46 } 47 48 // Test unit uniform 49 void test_unit_uniform(uint64_t* seed) 50 { 51 int n = 1000 * 1000; 52 double* unit_uniform_array = malloc(sizeof(double) * n); 53 54 for (int i = 0; i < n; i++) { 55 unit_uniform_array[i] = sample_unit_uniform(seed); 56 } 57 58 struct array_expectations expectations = { 59 .array = unit_uniform_array, 60 .n = n, 61 .name = "unit uniform", 62 .expected_mean = 0.5, 63 .expected_std = sqrt(1.0 / 12.0), 64 .tolerance = TOLERANCE, 65 }; 66 67 test_array_expectations(expectations); 68 free(unit_uniform_array); 69 } 70 71 // Test uniforms 72 void test_uniform(double start, double end, uint64_t* seed) 73 { 74 int n = 1000 * 1000; 75 double* uniform_array = malloc(sizeof(double) * n); 76 77 for (int i = 0; i < n; i++) { 78 uniform_array[i] = sample_uniform(start, end, seed); 79 } 80 81 char* name = malloc(MAX_NAME_LENGTH * sizeof(char)); 82 snprintf(name, MAX_NAME_LENGTH, "[%f, %f] uniform", start, end); 83 struct array_expectations expectations = { 84 .array = uniform_array, 85 .n = n, 86 .name = name, 87 .expected_mean = (start + end) / 2, 88 .expected_std = sqrt(1.0 / 12.0) * fabs(end - start), 89 .tolerance = fabs(end - start) * TOLERANCE, 90 }; 91 92 test_array_expectations(expectations); 93 free(name); 94 free(uniform_array); 95 } 96 97 // Test unit normal 98 void test_unit_normal(uint64_t* seed) 99 { 100 int n = 1000 * 1000; 101 double* unit_normal_array = malloc(sizeof(double) * n); 102 103 for (int i = 0; i < n; i++) { 104 unit_normal_array[i] = sample_unit_normal(seed); 105 } 106 107 struct array_expectations expectations = { 108 .array = unit_normal_array, 109 .n = n, 110 .name = "unit normal", 111 .expected_mean = 0, 112 .expected_std = 1, 113 .tolerance = TOLERANCE, 114 }; 115 116 test_array_expectations(expectations); 117 free(unit_normal_array); 118 } 119 120 // Test normal 121 void test_normal(double mean, double std, uint64_t* seed) 122 { 123 int n = 10 * 1000 * 1000; 124 double* normal_array = malloc(sizeof(double) * n); 125 126 for (int i = 0; i < n; i++) { 127 normal_array[i] = sample_normal(mean, std, seed); 128 } 129 130 char* name = malloc(MAX_NAME_LENGTH * sizeof(char)); 131 snprintf(name, MAX_NAME_LENGTH, "normal(%f, %f)", mean, std); 132 struct array_expectations expectations = { 133 .array = normal_array, 134 .n = n, 135 .name = name, 136 .expected_mean = mean, 137 .expected_std = std, 138 .tolerance = TOLERANCE, 139 }; 140 141 test_array_expectations(expectations); 142 free(name); 143 free(normal_array); 144 } 145 146 // Test lognormal 147 void test_lognormal(double logmean, double logstd, uint64_t* seed) 148 { 149 int n = 10 * 1000 * 1000; 150 double* lognormal_array = malloc(sizeof(double) * n); 151 152 for (int i = 0; i < n; i++) { 153 lognormal_array[i] = sample_lognormal(logmean, logstd, seed); 154 } 155 156 char* name = malloc(MAX_NAME_LENGTH * sizeof(char)); 157 snprintf(name, MAX_NAME_LENGTH, "lognormal(%f, %f)", logmean, logstd); 158 struct array_expectations expectations = { 159 .array = lognormal_array, 160 .n = n, 161 .name = name, 162 .expected_mean = exp(logmean + pow(logstd, 2) / 2), 163 .expected_std = sqrt((exp(pow(logstd, 2)) - 1) * exp(2 * logmean + pow(logstd, 2))), 164 .tolerance = TOLERANCE, 165 }; 166 167 test_array_expectations(expectations); 168 free(name); 169 free(lognormal_array); 170 } 171 172 // Test lognormal to 173 void test_to(double low, double high, uint64_t* seed) 174 { 175 int n = 10 * 1000 * 1000; 176 double* lognormal_array = malloc(sizeof(double) * n); 177 178 for (int i = 0; i < n; i++) { 179 lognormal_array[i] = sample_to(low, high, seed); 180 } 181 182 183 char* name = malloc(MAX_NAME_LENGTH * sizeof(char)); 184 snprintf(name, MAX_NAME_LENGTH, "to(%f, %f)", low, high); 185 186 const double NORMAL95CONFIDENCE = 1.6448536269514722; 187 double loglow = logf(low); 188 double loghigh = logf(high); 189 double logmean = (loglow + loghigh) / 2; 190 double logstd = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE); 191 192 struct array_expectations expectations = { 193 .array = lognormal_array, 194 .n = n, 195 .name = name, 196 .expected_mean = exp(logmean + pow(logstd, 2) / 2), 197 .expected_std = sqrt((exp(pow(logstd, 2)) - 1) * exp(2 * logmean + pow(logstd, 2))), 198 .tolerance = TOLERANCE, 199 }; 200 201 test_array_expectations(expectations); 202 free(name); 203 free(lognormal_array); 204 } 205 206 // Test beta 207 208 void test_beta(double a, double b, uint64_t* seed) 209 { 210 int n = 10 * 1000 * 1000; 211 double* beta_array = malloc(sizeof(double) * n); 212 213 for (int i = 0; i < n; i++) { 214 beta_array[i] = sample_beta(a, b, seed); 215 } 216 217 char* name = malloc(MAX_NAME_LENGTH * sizeof(char)); 218 snprintf(name, MAX_NAME_LENGTH, "beta(%f, %f)", a, b); 219 struct array_expectations expectations = { 220 .array = beta_array, 221 .n = n, 222 .name = name, 223 .expected_mean = a / (a + b), 224 .expected_std = sqrt((a * b) / (pow(a + b, 2) * (a + b + 1))), 225 .tolerance = TOLERANCE, 226 }; 227 228 test_array_expectations(expectations); 229 free(name); 230 } 231 232 int main() 233 { 234 // set randomness seed 235 uint64_t* seed = malloc(sizeof(uint64_t)); 236 *seed = 1000; // xorshift can't start with a seed of 0 237 238 printf("Testing unit uniform\n"); 239 test_unit_uniform(seed); 240 241 printf("Testing small uniforms\n"); 242 for (int i = 0; i < 100; i++) { 243 double start = sample_uniform(-10, 10, seed); 244 double end = sample_uniform(-10, 10, seed); 245 if (end > start) { 246 test_uniform(start, end, seed); 247 } 248 } 249 250 printf("Testing wide uniforms\n"); 251 for (int i = 0; i < 100; i++) { 252 double start = sample_uniform(-1000 * 1000, 1000 * 1000, seed); 253 double end = sample_uniform(-1000 * 1000, 1000 * 1000, seed); 254 if (end > start) { 255 test_uniform(start, end, seed); 256 } 257 } 258 259 printf("Testing unit normal\n"); 260 test_unit_normal(seed); 261 262 printf("Testing small normals\n"); 263 for (int i = 0; i < 100; i++) { 264 double mean = sample_uniform(-10, 10, seed); 265 double std = sample_uniform(0, 10, seed); 266 if (std > 0) { 267 test_normal(mean, std, seed); 268 } 269 } 270 271 printf("Testing larger normals\n"); 272 for (int i = 0; i < 100; i++) { 273 double mean = sample_uniform(-1000 * 1000, 1000 * 1000, seed); 274 double std = sample_uniform(0, 1000 * 1000, seed); 275 if (std > 0) { 276 test_normal(mean, std, seed); 277 } 278 } 279 280 printf("Testing smaller lognormals\n"); 281 for (int i = 0; i < 100; i++) { 282 double mean = sample_uniform(-1, 1, seed); 283 double std = sample_uniform(0, 1, seed); 284 if (std > 0) { 285 test_lognormal(mean, std, seed); 286 } 287 } 288 289 printf("Testing larger lognormals\n"); 290 for (int i = 0; i < 100; i++) { 291 double mean = sample_uniform(-1, 5, seed); 292 double std = sample_uniform(0, 5, seed); 293 if (std > 0) { 294 test_lognormal(mean, std, seed); 295 } 296 } 297 298 printf("Testing lognormals — sample_to(low, high) syntax\n"); 299 for (int i = 0; i < 100; i++) { 300 double low = sample_uniform(0, 1000 * 1000, seed); 301 double high = sample_uniform(0, 1000 * 1000, seed); 302 if (low < high) { 303 test_to(low, high, seed); 304 } 305 } 306 // Bonus example 307 test_to(10, 10 * 1000, seed); 308 309 printf("Testing beta distribution\n"); 310 for (int i = 0; i < 100; i++) { 311 double a = sample_uniform(0, 1000, seed); 312 double b = sample_uniform(0, 1000, seed); 313 if ((a > 0) && (b > 0)) { 314 test_beta(a, b, seed); 315 } 316 } 317 318 printf("Testing larger beta distributions\n"); 319 for (int i = 0; i < 100; i++) { 320 double a = sample_uniform(0, 1000 * 1000, seed); 321 double b = sample_uniform(0, 1000 * 1000, seed); 322 if ((a > 0) && (b > 0)) { 323 test_beta(a, b, seed); 324 } 325 } 326 327 free(seed); 328 }