| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright © 2025, Niklas Haas | ||
| 3 | * All rights reserved. | ||
| 4 | * | ||
| 5 | * Redistribution and use in source and binary forms, with or without | ||
| 6 | * modification, are permitted provided that the following conditions are met: | ||
| 7 | * | ||
| 8 | * 1. Redistributions of source code must retain the above copyright notice, this | ||
| 9 | * list of conditions and the following disclaimer. | ||
| 10 | * | ||
| 11 | * 2. Redistributions in binary form must reproduce the above copyright notice, | ||
| 12 | * this list of conditions and the following disclaimer in the documentation | ||
| 13 | * and/or other materials provided with the distribution. | ||
| 14 | * | ||
| 15 | * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND | ||
| 16 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED | ||
| 17 | * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE | ||
| 18 | * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR | ||
| 19 | * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES | ||
| 20 | * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; | ||
| 21 | * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND | ||
| 22 | * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT | ||
| 23 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | ||
| 24 | * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | ||
| 25 | */ | ||
| 26 | |||
| 27 | #include <math.h> | ||
| 28 | #include <stdlib.h> | ||
| 29 | |||
| 30 | #include "stats.h" | ||
| 31 | |||
| 32 | ✗ | CheckasmVar checkasm_var_scale(CheckasmVar a, double s) | |
| 33 | { | ||
| 34 | /* = checkasm_var_mul(a, checkasm_var_const(b)) */ | ||
| 35 | ✗ | return (CheckasmVar) { | |
| 36 | ✗ | .lmean = a.lmean + log(s), | |
| 37 | ✗ | .lvar = a.lvar, | |
| 38 | }; | ||
| 39 | } | ||
| 40 | |||
| 41 | ✗ | CheckasmVar checkasm_var_pow(CheckasmVar a, double exp) | |
| 42 | { | ||
| 43 | ✗ | return (CheckasmVar) { | |
| 44 | ✗ | .lmean = a.lmean * exp, | |
| 45 | ✗ | .lvar = a.lvar * exp * exp, | |
| 46 | }; | ||
| 47 | } | ||
| 48 | |||
| 49 | ✗ | CheckasmVar checkasm_var_add(const CheckasmVar a, const CheckasmVar b) | |
| 50 | { | ||
| 51 | /* Approximation assuming independent log-normal distributions */ | ||
| 52 | ✗ | const double ma = exp(a.lmean + 0.5 * a.lvar); | |
| 53 | ✗ | const double mb = exp(b.lmean + 0.5 * b.lvar); | |
| 54 | ✗ | const double va = (exp(a.lvar) - 1.0) * exp(2.0 * a.lmean + a.lvar); | |
| 55 | ✗ | const double vb = (exp(b.lvar) - 1.0) * exp(2.0 * b.lmean + b.lvar); | |
| 56 | ✗ | const double m = ma + mb; | |
| 57 | ✗ | const double v = va + vb; | |
| 58 | ✗ | return (CheckasmVar) { | |
| 59 | ✗ | .lmean = log(m * m / sqrt(v + m * m)), | |
| 60 | ✗ | .lvar = log(1.0 + v / (m * m)), | |
| 61 | }; | ||
| 62 | } | ||
| 63 | |||
| 64 | ✗ | CheckasmVar checkasm_var_sub(CheckasmVar a, CheckasmVar b) | |
| 65 | { | ||
| 66 | ✗ | const double ma = exp(a.lmean + 0.5 * a.lvar); | |
| 67 | ✗ | const double mb = exp(b.lmean + 0.5 * b.lvar); | |
| 68 | ✗ | const double va = (exp(a.lvar) - 1.0) * exp(2.0 * a.lmean + a.lvar); | |
| 69 | ✗ | const double vb = (exp(b.lvar) - 1.0) * exp(2.0 * b.lmean + b.lvar); | |
| 70 | ✗ | const double m = fmax(ma - mb, 1e-30); /* avoid negative mean */ | |
| 71 | ✗ | const double v = va + vb; | |
| 72 | ✗ | return (CheckasmVar) { | |
| 73 | ✗ | .lmean = log(m * m / sqrt(v + m * m)), | |
| 74 | ✗ | .lvar = log(1.0 + v / (m * m)), | |
| 75 | }; | ||
| 76 | } | ||
| 77 | |||
| 78 | ✗ | CheckasmVar checkasm_var_mul(CheckasmVar a, CheckasmVar b) | |
| 79 | { | ||
| 80 | ✗ | return (CheckasmVar) { | |
| 81 | ✗ | .lmean = a.lmean + b.lmean, | |
| 82 | ✗ | .lvar = a.lvar + b.lvar, | |
| 83 | }; | ||
| 84 | } | ||
| 85 | |||
| 86 | ✗ | CheckasmVar checkasm_var_inv(CheckasmVar a) | |
| 87 | { | ||
| 88 | ✗ | return (CheckasmVar) { | |
| 89 | ✗ | .lmean = -a.lmean, | |
| 90 | ✗ | .lvar = a.lvar, | |
| 91 | }; | ||
| 92 | } | ||
| 93 | |||
| 94 | ✗ | CheckasmVar checkasm_var_div(CheckasmVar a, CheckasmVar b) | |
| 95 | { | ||
| 96 | ✗ | return (CheckasmVar) { | |
| 97 | ✗ | .lmean = a.lmean - b.lmean, | |
| 98 | ✗ | .lvar = a.lvar + b.lvar, | |
| 99 | }; | ||
| 100 | } | ||
| 101 | |||
| 102 | ✗ | CheckasmVar checkasm_stats_estimate(const CheckasmStats *const stats) | |
| 103 | { | ||
| 104 | ✗ | if (!stats->nb_samples) | |
| 105 | ✗ | return checkasm_var_const(0.0); | |
| 106 | |||
| 107 | /* Compute mean and variance */ | ||
| 108 | ✗ | double sum = 0.0, sum2 = 0.0, sum_w2 = 0.0; | |
| 109 | ✗ | int count = 0; | |
| 110 | ✗ | for (int i = 0; i < stats->nb_samples; i++) { | |
| 111 | ✗ | const CheckasmSample s = stats->samples[i]; | |
| 112 | ✗ | const double x = log((double) s.sum) - log((double) s.count); | |
| 113 | ✗ | sum += x * s.count; | |
| 114 | ✗ | sum2 += x * x * s.count; | |
| 115 | ✗ | sum_w2 += (double) s.count * s.count; | |
| 116 | ✗ | count += s.count; | |
| 117 | } | ||
| 118 | |||
| 119 | ✗ | assert(count > 0); | |
| 120 | ✗ | const double mean = sum / count; | |
| 121 | ✗ | const double denom = count - sum_w2 / count; | |
| 122 | double var; | ||
| 123 | ✗ | if (denom > 0.0) { | |
| 124 | ✗ | var = fmax(sum2 - count * mean * mean, 0.0) / denom; | |
| 125 | } else { | ||
| 126 | /* Lower bound on the variance predicted by the sample count alone */ | ||
| 127 | ✗ | var = 1.0 / count; | |
| 128 | } | ||
| 129 | |||
| 130 | ✗ | return (CheckasmVar) { .lmean = mean, .lvar = var }; | |
| 131 | } | ||
| 132 |