| line | % | coverage | branch |
|---|
| 30 | 50 | T | F | * PERL_NO_GET_CONTEXT is therefore not a concern here. |
| 40 | 100 | T | F | z = (z ^ (z >> 27)) * UINT64_C(0x94d049bb133111eb); |
| 52 | 100 | T | F | s = (uint64_t)time(NULL) ^ ((uint64_t)getpid() << 32); |
| 73 | 100 | T | F | (void)seedDrand01((Rand_seed_t)Perl_seed(aTHX)); \ |
| 79 | 100 | T | F | // Helpers for Random Number Generation |
| 81 | 100 | T | F | #ifndef M_PI |
| 50 | T | F | #ifndef M_PI |
| 83 | 100 | T | F | #endif |
| 87 | 100 | T | F | if (df <= 0.0) return 0.0; |
| 93 | 100 | T | F | double root_half = 0.70710678118654752440; // 1 / sqrt(2) |
| 94 | 100 | T | F | |
| 97 | 100 | T | F | double w = u / (1.0 - u); |
| 98 | 100 | T | F | // Scaled Chi-distribution log-density |
| 111 | 100 | T | F | // Ranking helper with tie adjustment (matches R's tie handling) |
| 113 | 100 | T | F | static int compare_rank(const void *restrict a, const void *restrict b) { |
| 117 | 100 | T | F | |
| 124 | 100 | T | F | for (size_t i = 0; i < n; i++) { |
| 126 | 100 | T | F | items[i].idx = i; |
| 127 | 100 | T | F | } |
| 135 | 100 | T | F | i = j; |
| 143 | 100 | T | F | static size_t generate_binomial(const size_t size, const double prob) { |
| 148 | 50 | T | F | for (size_t i = 0; i < size; i++) { |
| 0 | T | F | for (size_t i = 0; i < size; i++) { |
| 149 | 50 | T | F | if (Drand01() <= prob) successes++; |
| 150 | 100 | T | F | } |
| 153 | 50 | T | F | // Helper: log combination |
| 156 | 100 | T | F | } |
| 158 | 100 | T | F | // Log-space tails for non-central hypergeometric |
| 161 | 100 | T | F | |
| 168 | 100 | T | F | for(size_t k = 0; k <= max_x - min_x; ++k) { |
| 170 | 100 | T | F | } |
| 179 | 100 | T | F | } |
| 180 | 100 | T | F | } |
| 181 | 50 | T | F | |
| 183 | 50 | T | F | static void calculate_exact_stats(size_t a, size_t b, size_t c, size_t d, double conf_level, const char*restrict alt, double *restrict mle_or, double *restrict ci_low, double *restrict ci_high) { |
| 188 | 100 | T | F | |
| 189 | 100 | T | F | bool is_less = (strcmp(alt, "less") == 0); |
| 191 | 100 | T | F | |
| 198 | 100 | T | F | // MLE |
| 199 | 100 | T | F | if (a == min_x && a == max_x) *mle_or = 1.0; |
| 200 | 100 | T | F | else if (a == min_x) *mle_or = 0.0; |
| 202 | 50 | T | F | else { |
| 207 | 100 | T | F | for(size_t k = 0; k <= max_x - min_x; ++k) { |
| 208 | 100 | T | F | double d_val = logdc[k] + log_mid * (min_x + k); |
| 210 | 100 | T | F | } |
| 221 | 100 | T | F | } |
| 226 | 100 | T | F | *ci_high = INFINITY; |
| 232 | 100 | T | F | double log_low = -100.0, log_high = 100.0, best = 1.0, best_err = 1e9, lt, ut; |
| 233 | 100 | T | F | for (unsigned short int i = 0; i < 1000; i++) { |
| 234 | 100 | T | F | double log_mid = 0.5 * (log_low + log_high); |
| 235 | 100 | T | F | double mid = exp(log_mid); |
| 239 | 100 | T | F | if (ut > target_alpha) log_high = log_mid; |
| 241 | 100 | T | F | if (log_high - log_low < 1e-15) break; |
| 246 | 50 | T | F | |
| 261 | 100 | T | F | } |
| 266 | 100 | T | F | } |
| 269 | 100 | T | F | static double exact_p_value(size_t a, size_t b, size_t c, size_t d, const char* alt) { |
| 50 | T | F | static double exact_p_value(size_t a, size_t b, size_t c, size_t d, const char* alt) { |
| 272 | 100 | T | F | size_t max_x = (r1 < c1) ? r1 : c1; |
| 281 | 100 | T | F | |
| 282 | 100 | T | F | if (strcmp(alt, "less") == 0) { |
| 283 | 100 | T | F | for(size_t x = min_x; x <= a; ++x) p_val += exp(logdc[x - min_x]); |
| 100 | T | F | for(size_t x = min_x; x <= a; ++x) p_val += exp(logdc[x - min_x]); |
| 286 | 100 | T | F | } else { |
| 299 | 100 | T | F | * Helpers for lm Linear Regression: OLS Matrix Math & Formula Parsing |
| 301 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 305 | 50 | T | F | */ |
| 307 | 50 | T | F | int rank = 0; |
| 50 | T | F | int rank = 0; |
| 50 | T | F | int rank = 0; |
| 312 | 50 | T | F | aliased[k] = FALSE; |
| 100 | T | F | aliased[k] = FALSE; |
| 313 | 100 | T | F | orig_diag[k] = A[k * n + k]; |
| 322 | 50 | T | F | for (size_t i = 0; i < n; i++) { |
| 325 | 100 | T | F | } |
| 328 | 0 | T | F | rank++; |
| 0 | T | F | rank++; |
| 0 | T | F | rank++; |
| 331 | 0 | T | F | for (size_t j = 0; j < n; j++) A[k * n + j] *= pivot; |
| 340 | 50 | T | F | } |
| 50 | T | F | } |
| 344 | 100 | T | F | } |
| 350 | 50 | T | F | val = hv_fetch(row_hashes[i], var, strlen(var), 0); |
| 50 | T | F | val = hv_fetch(row_hashes[i], var, strlen(var), 0); |
| 353 | 50 | T | F | val = av_fetch(av, 0, 0); |
| 355 | 0 | T | F | } else if (data_hoa) { |
| 359 | 0 | T | F | val = av_fetch(av, i, 0); |
| 366 | 0 | T | F | return NAN; // Catch undef/missing keys |
| 367 | 0 | T | F | } |
| 376 | 100 | T | F | av_push(cols, newSVsv(hv_iterkeysv(entry))); |
| 378 | 100 | T | F | } else if (row_hashes && n > 0 && row_hashes[0]) { |
| 380 | 100 | T | F | HE *restrict entry; |
| 50 | T | F | HE *restrict entry; |
| 50 | T | F | HE *restrict entry; |
| 384 | 50 | T | F | } |
| 386 | 50 | T | F | } |
| 50 | T | F | } |
| 50 | T | F | } |
| 391 | 100 | T | F | |
| 50 | T | F | |
| 392 | 100 | T | F | char *restrict term_cpy = savepv(term); |
| 402 | 50 | T | F | } |
| 404 | 0 | T | F | char *restrict end = strrchr(term_cpy, ')'); |
| 0 | T | F | char *restrict end = strrchr(term_cpy, ')'); |
| 0 | T | F | char *restrict end = strrchr(term_cpy, ')'); |
| 408 | 50 | T | F | int power = 1; |
| 410 | 50 | T | F | *caret = '\0'; |
| 50 | T | F | *caret = '\0'; |
| 50 | T | F | *caret = '\0'; |
| 415 | 50 | T | F | |
| 50 | T | F | |
| 430 | 100 | T | F | if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVAV) { |
| 431 | 50 | T | F | AV*restrict av = (AV*)SvRV(*val); |
| 447 | 100 | T | F | } |
| 448 | 100 | T | F | |
| 456 | 50 | T | F | val = av_fetch(av, 0, 0); |
| 457 | 100 | T | F | } |
| 461 | 100 | T | F | AV*restrict av = (AV*)SvRV(*col); |
| 464 | 100 | T | F | } |
| 100 | T | F | } |
| 467 | 100 | T | F | } |
| 477 | 100 | T | F | // Comparator for qsort |
| 483 | 50 | T | F | return ((PVal*)a)->orig_idx - ((PVal*)b)->orig_idx; |
| 497 | 100 | T | F | if (diff < 0) return -1; |
| 498 | 100 | T | F | if (diff > 0) return 1; |
| 501 | 100 | T | F | |
| 50 | T | F | |
| 502 | 100 | T | F | /* Compute 1-based average ranks with tie-breaking into out[]. |
| 503 | 50 | T | F | * in[] is not modified. */ |
| 504 | 50 | T | F | static void rank_data(const double *restrict in, double *restrict out, size_t n) { |
| 509 | 50 | T | F | |
| 517 | 100 | T | F | for (size_t k = i; k <= j; k++) out[ri[k].idx] = avg; |
| 519 | 50 | T | F | } |
| 50 | T | F | } |
| 526 | 100 | T | F | double sx = 0, sy = 0, sxy = 0, sx2 = 0, sy2 = 0; |
| 542 | 50 | T | F | * x, T_y = pairs tied only on y. Joint ties (both zero) are excluded |
| 544 | 50 | T | F | * Returns NAN when the denominator is zero. */ |
| 548 | 50 | T | F | for (size_t j = i + 1; j < n; j++) { |
| 550 | 50 | T | F | int sy = (y[i] > y[j]) - (y[i] < y[j]); |
| 554 | 50 | T | F | else if (sx == sy) C++; |
| 556 | 50 | T | F | } |
| 558 | 100 | T | F | double denom = sqrt((double)(C + D + tie_x) * (double)(C + D + tie_y)); |
| 564 | 50 | T | F | * Allocates and frees temporary rank arrays internally for Spearman. */ |
| 565 | 100 | T | F | static double compute_cor(const double *restrict x, const double *restrict y, |
| 567 | 100 | T | F | if (strcmp(method, "spearman") == 0) { |
| 574 | 100 | T | F | return r; |
| 100 | T | F | return r; |
| 575 | 100 | T | F | } |
| 50 | T | F | } |
| 583 | 100 | T | F | #define MAX_ITER 500 |
| 586 | 50 | T | F | |
| 589 | 50 | T | F | double aa, c, d, del, h, qab, qam, qap; |
| 592 | 100 | T | F | if (fabs(d) < FPMIN) d = FPMIN; |
| 597 | 100 | T | F | d = 1.0 + aa * d; |
| 609 | 0 | T | F | } |
| 618 | 50 | T | F | return 1.0 - bt * _incbeta_cf(b, a, 1.0 - x) / b; |
| 620 | 100 | T | F | |
| 625 | 100 | T | F | if (strcmp(alt, "greater") == 0) return (t > 0) ? 0.5 * prob_2tail : 1.0 - 0.5 * prob_2tail; |
| 626 | 100 | T | F | return prob_2tail; |
| 629 | 50 | T | F | // Bisection algorithm to find the inverse t-distribution (Critical t-value) |
| 50 | T | F | // Bisection algorithm to find the inverse t-distribution (Critical t-value) |
| 50 | T | F | // Bisection algorithm to find the inverse t-distribution (Critical t-value) |
| 633 | 100 | T | F | while (get_t_pvalue(high, df, "greater") > p_tail) { |
| 639 | 100 | T | F | for (unsigned short int i = 0; i < 100; i++) { |
| 100 | T | F | for (unsigned short int i = 0; i < 100; i++) { |
| 643 | 100 | T | F | low = mid; |
| 50 | T | F | low = mid; |
| 648 | 50 | T | F | } |
| 653 | 100 | T | F | double da = *(const double*restrict)a; |
| 655 | 100 | T | F | return (da > db) - (da < db); |
| 658 | 50 | T | F | static size_t calculate_sturges_bins(size_t n) { |
| 681 | 100 | T | F | size_t idx = (size_t)((val - min_val) / step); |
| 687 | 100 | T | F | /* If value is exactly on or slightly below the lower boundary of the assigned bin, |
| 691 | 100 | T | F | } |
| 709 | 100 | T | F | } |
| 724 | 100 | T | F | double a[4] = {2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637}; |
| 50 | T | F | double a[4] = {2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637}; |
| 50 | T | F | double a[4] = {2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637}; |
| 727 | 100 | T | F | 0.0276438810333863, 0.0038405729373609, 0.0003951896511919, |
| 728 | 100 | T | F | 0.0000321767881768, 0.0000002888167364, 0.0000003960315187}; |
| 730 | 100 | T | F | y = p - 0.5; |
| 735 | 100 | T | F | } else { |
| 100 | T | F | } else { |
| 100 | T | F | } else { |
| 751 | 50 | T | F | * i.e. rho ≥ rho_obs) and how many yield S ≥ s_obs ("upper tail"). |
| 752 | 50 | T | F | * |
| 754 | 50 | T | F | * Valid up to n = 9 (362 880 iterations — negligible cost). |
| 755 | 50 | T | F | * ----------------------------------------------------------------------- */ |
| 764 | 100 | T | F | double s_ = 0.0; \ |
| 767 | 100 | T | F | s_ += d_ * d_; \ |
| 769 | 100 | T | F | if (s_ <= s_obs + 1e-9) count_le++; \ |
| 771 | 100 | T | F | total++; \ |
| 773 | 100 | T | F | |
| 100 | T | F | |
| 784 | 50 | T | F | } |
| 785 | 50 | T | F | TALLY_PERM(); |
| 787 | 100 | T | F | k = 1; |
| 789 | 100 | T | F | c[k] = 0; |
| 791 | 50 | T | F | } |
| 792 | 100 | T | F | } |
| 794 | 50 | T | F | |
| 795 | 50 | T | F | Safefree(perm); Safefree(c); |
| 799 | 50 | T | F | double p_ge = (double)count_ge / (double)total; |
| 808 | 100 | T | F | * Exact Kendall p-value via Mahonian Numbers (Inversions distribution) |
| 810 | 50 | T | F | * ----------------------------------------------------------------------- */ |
| 816 | 100 | T | F | /* Build the distribution of inversions via DP */ |
| 817 | 100 | T | F | for (size_t i = 2; i <= n; i++) { |
| 819 | 100 | T | F | for (long k = 0; k <= max_inv; k++) next_dp[k] = 0.0; |
| 825 | 100 | T | F | } |
| 830 | 100 | T | F | dp = next_dp; |
| 834 | 100 | T | F | if (i_obs < 0) i_obs = 0; |
| 836 | 100 | T | F | double p_le = 0.0; /* P(S <= S_obs) */ |
| 839 | 100 | T | F | for (long k = 0; k <= i_obs; k++) p_ge += dp[k]; |
| 844 | 100 | T | F | double p = 2.0 * (p_ge < p_le ? p_ge : p_le); |
| 847 | 100 | T | F | // F-distribution Cumulative Distribution Function P(F <= f) |
| 863 | 50 | T | F | } |
| 50 | T | F | } |
| 866 | 100 | T | F | for (size_t i = r; i < n; i++) { |
| 867 | 50 | T | F | if (fabs(X[i][k]) > max_val) max_val = fabs(X[i][k]); |
| 874 | 50 | T | F | double norm = 0; |
| 876 | 100 | T | F | X[i][k] /= max_val; |
| 100 | T | F | X[i][k] /= max_val; |
| 50 | T | F | X[i][k] /= max_val; |
| 50 | T | F | X[i][k] /= max_val; |
| 880 | 100 | T | F | double s = (X[r][k] > 0) ? -norm : norm; |
| 882 | 100 | T | F | X[r][k] = s * max_val; |
| 883 | 100 | T | F | |
| 899 | 100 | T | F | rank_map[k] = r; // Map original column index to orthogonal row index |
| 900 | 100 | T | F | r++; |
| 901 | 50 | T | F | } |
| 912 | 50 | T | F | static bool contains_nondigit(SV *restrict sv) { |
| 50 | T | F | static bool contains_nondigit(SV *restrict sv) { |
| 913 | 50 | T | F | if (!sv || !SvOK(sv)) return 0; |
| 916 | 100 | T | F | for (size_t i = 0; i < len; i++) { |
| 920 | 100 | T | F | } |
| 933 | 50 | T | F | if (*p == '"') { |
| 937 | 50 | T | F | PerlIO_putc(fh, *p); |
| 939 | 50 | T | F | } |
| 943 | 100 | T | F | } |
| 951 | 50 | T | F | if (row[i]) { |
| 952 | 50 | T | F | print_str_quoted(fh, row[i], sep); |
| 963 | 50 | T | F | if (x == 0.0) return 1.0; |
| 50 | T | F | if (x == 0.0) return 1.0; |
| 964 | 50 | T | F | |
| 966 | 100 | T | F | if (x < a + 1.0) { |
| 977 | 100 | T | F | |
| 978 | 50 | T | F | // Continued fraction for x >= a + 1 |
| 983 | 100 | T | F | while (i < 10000) { // Safety bound |
| 984 | 100 | T | F | double an = -i * (i - a); |
| 985 | 100 | T | F | b += 2.0; |
| 989 | 100 | T | F | if (fabs(c) < 1e-30) c = 1e-30; |
| 1003 | 50 | T | F | return igamc((double)df / 2.0, stat / 2.0); |
| 1004 | 100 | T | F | } |
| 1009 | 100 | T | F | #endif |
| 1010 | 100 | T | F | |
| 1014 | 100 | T | F | if (k > n / 2) k = n - k; |
| 1030 | 50 | T | F | long double *restrict w = (long double *)safecalloc(max_u + 1, sizeof(long double)); |
| 1035 | 100 | T | F | for (int i = max_u; i >= j + m; i--) w[i] -= w[i - j - m]; |
| 1037 | 100 | T | F | |
| 100 | T | F | |
| 1039 | 100 | T | F | for (int i = 0; i <= k; i++) cum_p += w[i]; |
| 1041 | 100 | T | F | long double total = choose_comb(m + n, n); |
| 1059 | 50 | T | F | for (int i = 1; i <= n; i++) { |
| 1060 | 50 | T | F | for (int j = max_v; j >= i; j--) w[j] += w[j - i]; |
| 1062 | 100 | T | F | |
| 1070 | 0 | T | F | return result; |
| 1071 | 0 | T | F | } |
| 1073 | 0 | T | F | static int cmp_rank_info(const void *a, const void *b) { |
| 1078 | 0 | T | F | |
| 1082 | 0 | T | F | size_t i = 0; |
| 1087 | 0 | T | F | while (j < n && ri[j].val == ri[i].val) j++; |
| 1092 | 0 | T | F | i = j; |
| 1105 | 100 | T | F | #endif |
| 1106 | 100 | T | F | |
| 1108 | 100 | T | F | static double r_pow_di(double x, int n) { |
| 1115 | 100 | T | F | |
| 1116 | 100 | T | F | // Two-sample two-sided asymptotic distribution |
| 1124 | 100 | T | F | int k_max = (int) sqrt(2.0 - log(tol)); |
| 1125 | 100 | T | F | double w = log(x); |
| 1131 | 50 | T | F | p = s / M_1_SQRT_2PI; |
| 1132 | 0 | T | F | if(!lower) p = 1.0 - p; |
| 1146 | 100 | T | F | k++; |
| 1147 | 100 | T | F | } |
| 1148 | 100 | T | F | p = new_val; |
| 1152 | 100 | T | F | |
| 1156 | 50 | T | F | for(unsigned int j = 0; j < m; j++) { |
| 1158 | 100 | T | F | for(unsigned int k = 0; k < m; k++) s += A[i * m + k] * B[k * m + j]; |
| 1159 | 100 | T | F | C[i * m + j] = s; |
| 1160 | 100 | T | F | } |
| 1161 | 100 | T | F | } |
| 1170 | 100 | T | F | m_power(A, eA, V, eV, m, n / 2); |
| 1172 | 50 | T | F | m_multiply(V, V, B, m); |
| 1191 | 100 | T | F | int m = 2 * k - 1; |
| 50 | T | F | int m = 2 * k - 1; |
| 1193 | 50 | T | F | double *restrict H = (double*) safecalloc(m * m, sizeof(double)); |
| 100 | T | F | double *restrict H = (double*) safecalloc(m * m, sizeof(double)); |
| 100 | T | F | double *restrict H = (double*) safecalloc(m * m, sizeof(double)); |
| 1194 | 50 | T | F | double *restrict Q = (double*) safecalloc(m * m, sizeof(double)); |
| 1197 | 100 | T | F | for(int j = 0; j < m; j++) { |
| 100 | T | F | for(int j = 0; j < m; j++) { |
| 1198 | 100 | T | F | if(i - j + 1 < 0) H[i * m + j] = 0; |
| 100 | T | F | if(i - j + 1 < 0) H[i * m + j] = 0; |
| 1204 | 100 | T | F | H[(m - 1) * m + i] -= r_pow_di(h, (m - i)); |
| 1205 | 100 | T | F | } |
| 1206 | 100 | T | F | H[(m - 1) * m] += ((2 * h - 1 > 0) ? r_pow_di(2 * h - 1, m) : 0); |
| 1215 | 100 | T | F | |
| 1225 | 100 | T | F | } |
| 1226 | 100 | T | F | } |
| 1229 | 100 | T | F | Safefree(Q); |
| 1230 | 100 | T | F | return s; |
| 1231 | 100 | T | F | } |
| 1232 | 100 | T | F | |
| 1247 | 50 | T | F | while(i < nx && x[i] <= val) i++; |
| 1255 | 50 | T | F | if (-diff > max_d_minus) max_d_minus = -diff; |
| 0 | T | F | if (-diff > max_d_minus) max_d_minus = -diff; |
| 1268 | 50 | T | F | |
| 1269 | 50 | T | F | // Evaluate the exact 2-sample probability |
| 1272 | 100 | T | F | double *restrict u = (double *) safecalloc(n + 1, sizeof(double)); |
| 1275 | 50 | T | F | for(unsigned int j = 1; j <= n; j++) { |
| 1279 | 50 | T | F | for(unsigned int i = 1; i <= m; i++) { |
| 1283 | 100 | T | F | else { |
| 1288 | 100 | T | F | } |
| 1304 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 1309 | 50 | T | F | // Default: 1 - P(T < qu) |
| 1310 | 100 | T | F | // Ensure exact_pnt is using a convergence tolerance of at least 1e-15 |
| 50 | T | F | // Ensure exact_pnt is using a convergence tolerance of at least 1e-15 |
| 1313 | 50 | T | F | } |
| 1320 | 100 | T | F | double low = 0.0, high = 1.0; |
| 1323 | 50 | T | F | low = high; |
| 1324 | 50 | T | F | high *= 2.0; |
| 1325 | 50 | T | F | if (high > 1e100) break; /* Fallback limit */ |
| 1326 | 0 | T | F | } |
| 1329 | 50 | T | F | for (unsigned short int i = 0; i < 150; i++) { |
| 1333 | 50 | T | F | if (p_mid < p) { |
| 50 | T | F | if (p_mid < p) { |
| 50 | T | F | if (p_mid < p) { |
| 1341 | 100 | T | F | } |
| 100 | T | F | } |
| 50 | T | F | } |
| 1347 | 50 | T | F | * my $res = oneway_test(\%groups, var_equal => 1); |
| 1352 | 100 | T | F | * ── Mode 2: formula – response ~ factor ─────────────────────────────────── |
| 1354 | 50 | T | F | * my $res = oneway_test(\%data, formula => "yield ~ ctrl"); |
| 50 | T | F | * my $res = oneway_test(\%data, formula => "yield ~ ctrl"); |
| 50 | T | F | * my $res = oneway_test(\%data, formula => "yield ~ ctrl"); |
| 1363 | 50 | T | F | * oneway.test(Value ~ Group, data = my_data) |
| 100 | T | F | * oneway.test(Value ~ Group, data = my_data) |
| 50 | T | F | * oneway.test(Value ~ Group, data = my_data) |
| 1369 | 100 | T | F | * Hash ref with keys: |
| 1371 | 50 | T | F | * num_df => numerator degrees of freedom (k − 1) |
| 50 | T | F | * num_df => numerator degrees of freedom (k − 1) |
| 50 | T | F | * num_df => numerator degrees of freedom (k − 1) |
| 1376 | 50 | T | F | * n => total observations |
| 50 | T | F | * n => total observations |
| 1385 | 100 | T | F | /* ----------------------------------------------------------------------- |
| 1386 | 100 | T | F | * C HELPERS (place above "--- XS SECTION ---") |
| 1391 | 50 | T | F | double statistic; |
| 1392 | 50 | T | F | double num_df; |
| 1398 | 100 | T | F | double ms_within; /* ss_within / denom_df */ |
| 1399 | 100 | T | F | int k; /* number of groups */ |
| 1403 | 100 | T | F | |
| 1404 | 50 | T | F | /* ── c_oneway_test ─────────────────────────────────────────────────────── |
| 1407 | 50 | T | F | * sizes – n_i for each group (length k) |
| 50 | T | F | * sizes – n_i for each group (length k) |
| 1411 | 50 | T | F | * Mirrors R's oneway.test() arithmetic exactly. |
| 1418 | 0 | T | F | int var_equal) |
| 1427 | 50 | T | F | |
| 50 | T | F | |
| 1429 | 50 | T | F | IV total_n = 0; |
| 1432 | 100 | T | F | n_i[g] = (double)ng; |
| 1440 | 50 | T | F | double ss = 0.0; |
| 1441 | 100 | T | F | for (size_t i = 0; i < ng; i++) { |
| 1442 | 100 | T | F | double d = data[offset + i] - mean; |
| 1443 | 50 | T | F | ss += d * d; |
| 1445 | 100 | T | F | v_i[g] = ss / (double)(ng - 1); /* ng >= 2 guaranteed by caller */ |
| 1446 | 50 | T | F | offset += ng; |
| 1448 | 50 | T | F | |
| 1449 | 50 | T | F | res.n = total_n; |
| 1451 | 50 | T | F | /* grand mean (simple average over all obs; used only by classic branch) */ |
| 1452 | 50 | T | F | double grand_mean = 0.0; |
| 1454 | 50 | T | F | grand_mean /= (double)total_n; |
| 1464 | 0 | T | F | double dm = m_i[g] - grand_mean; |
| 1476 | 50 | T | F | } else { |
| 1477 | 50 | T | F | /* ── Welch one-way (heteroscedastic) ───────────────────────────── * |
| 1498 | 50 | T | F | } |
| 100 | T | F | } |
| 50 | T | F | } |
| 1503 | 50 | T | F | num += w_i[g] * dm * dm; |
| 100 | T | F | num += w_i[g] * dm * dm; |
| 50 | T | F | num += w_i[g] * dm * dm; |
| 1508 | 50 | T | F | /* unweighted SS for the output table */ |
| 1512 | 100 | T | F | ssbg += n_i[g] * dm * dm; |
| 1515 | 100 | T | F | res.ss_between = ssbg; |
| 1516 | 100 | T | F | res.ss_within = sswg; |
| 1517 | 100 | T | F | res.ms_between = (df1 > 0.0) ? ssbg / df1 : 0.0; |
| 1518 | 50 | T | F | res.ms_within = (res.denom_df > 0.0) ? sswg / res.denom_df : 0.0; |
| 1519 | 100 | T | F | Safefree(w_i); |
| 1520 | 50 | T | F | } |
| 1521 | 0 | T | F | /* upper-tail p-value P(F ≥ statistic) */ |
| 1524 | 50 | T | F | return res; |
| 1528 | 100 | T | F | * |
| 50 | T | F | * |
| 50 | T | F | * |
| 1532 | 50 | T | F | * Caller must Safefree() both *lhs and *rhs on success. |
| 1536 | 100 | T | F | { |
| 50 | T | F | { |
| 50 | T | F | { |
| 1544 | 100 | T | F | if (l_end < l_start) return 0; /* empty LHS */ |
| 100 | T | F | if (l_end < l_start) return 0; /* empty LHS */ |
| 1547 | 100 | T | F | const char *restrict r_start = tilde + 1; |
| 1549 | 50 | T | F | const char *restrict r_end = r_start + strlen(r_start) - 1; |
| 50 | T | F | const char *restrict r_end = r_start + strlen(r_start) - 1; |
| 50 | T | F | const char *restrict r_end = r_start + strlen(r_start) - 1; |
| 1555 | 100 | T | F | |
| 1557 | 50 | T | F | *rhs = (char *)safemalloc(rlen + 1); |
| 50 | T | F | *rhs = (char *)safemalloc(rlen + 1); |
| 50 | T | F | *rhs = (char *)safemalloc(rlen + 1); |
| 1563 | 50 | T | F | /* ── build_groups_from_formula ─────────────────────────────────────────── |
| 1564 | 50 | T | F | * |
| 1569 | 100 | T | F | * slots for both; actual group count returned via *out_k) |
| 100 | T | F | * slots for both; actual group count returned via *out_k) |
| 1572 | 50 | T | F | * |
| 1573 | 50 | T | F | * Group identity is the string representation of each label element |
| 1574 | 50 | T | F | * (SvPV_nolen), so integer 0 and string "0" are the same group. |
| 50 | T | F | * (SvPV_nolen), so integer 0 and string "0" are the same group. |
| 100 | T | F | * (SvPV_nolen), so integer 0 and string "0" are the same group. |
| 1576 | 100 | T | F | * factor level ordering from stack(). |
| 50 | T | F | * factor level ordering from stack(). |
| 1580 | 100 | T | F | #define OWT_MAX_GROUPS 1024 /* sane ceiling; ANOVA with >1024 groups is absurd */ |
| 1585 | 100 | T | F | AV *restrict label_av, |
| 1586 | 50 | T | F | double *restrict out_flat, |
| 1588 | 0 | T | F | size_t *restrict out_k, |
| 1592 | 50 | T | F | { |
| 1598 | 50 | T | F | "formula: response length (%"IVdf") != factor length (%"IVdf")", |
| 1599 | 50 | T | F | n, nl); |
| 100 | T | F | n, nl); |
| 1600 | 0 | T | F | return 0; |
| 1601 | 0 | T | F | } |
| 1605 | 50 | T | F | } |
| 1606 | 50 | T | F | |
| 1611 | 100 | T | F | IV *restrict obs_group = (IV *)safemalloc((size_t)n * sizeof(IV)); |
| 50 | T | F | IV *restrict obs_group = (IV *)safemalloc((size_t)n * sizeof(IV)); |
| 100 | T | F | IV *restrict obs_group = (IV *)safemalloc((size_t)n * sizeof(IV)); |
| 1615 | 100 | T | F | SV **lsv = av_fetch(label_av, i, 0); |
| 1617 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 1620 | 100 | T | F | for (size_t g = 0; g < ngroups; g++) { |
| 1622 | 50 | T | F | } |
| 50 | T | F | } |
| 50 | T | F | } |
| 1625 | 50 | T | F | snprintf(errbuf, errbuf_len, |
| 1629 | 50 | T | F | return 0; |
| 1633 | 50 | T | F | group_names[ngroups] = (char *)safemalloc(lablen + 1); |
| 1638 | 100 | T | F | } |
| 1645 | 100 | T | F | return 0; |
| 1646 | 100 | T | F | } |
| 1648 | 50 | T | F | /* count per-group sizes */ |
| 1649 | 50 | T | F | memset(out_sizes, 0, ngroups * sizeof(size_t)); |
| 1650 | 50 | T | F | for (unsigned i = 0; i < n; i++) out_sizes[obs_group[i]]++; |
| 50 | T | F | for (unsigned i = 0; i < n; i++) out_sizes[obs_group[i]]++; |
| 1651 | 50 | T | F | |
| 50 | T | F | |
| 1655 | 50 | T | F | snprintf(errbuf, errbuf_len, |
| 50 | T | F | snprintf(errbuf, errbuf_len, |
| 1659 | 50 | T | F | Safefree(group_names); Safefree(obs_group); |
| 1664 | 50 | T | F | /* ── fill flat output array in group order ─────────────────────── * |
| 1665 | 50 | T | F | * We compute a running write-offset per group, then scatter. * |
| 1667 | 50 | T | F | size_t *restrict write_pos = (size_t *)safemalloc(ngroups * sizeof(size_t)); |
| 1671 | 0 | T | F | |
| 1676 | 0 | T | F | out_flat[write_pos[g]++] = val; |
| 1677 | 0 | T | F | } |
| 0 | T | F | } |
| 1678 | 0 | T | F | |
| 1679 | 0 | T | F | *out_k = ngroups; |
| 1683 | 0 | T | F | if (out_names) { |
| 1684 | 0 | T | F | *out_names = group_names; /* caller takes ownership */ |
| 1689 | 50 | T | F | return 1; |
| 1705 | 50 | T | F | char *lhs = NULL, *rhs = NULL; |
| 1708 | 50 | T | F | char ** gnames = NULL; |
| 100 | T | F | char ** gnames = NULL; |
| 50 | T | F | char ** gnames = NULL; |
| 1711 | 50 | T | F | IV total_n = 0; |
| 1719 | 100 | T | F | SV *restrict val = ST(ai + 1); |
| 50 | T | F | SV *restrict val = ST(ai + 1); |
| 100 | T | F | SV *restrict val = ST(ai + 1); |
| 1722 | 100 | T | F | else if (strEQ(key, "formula")) |
| 1725 | 100 | T | F | /* validate data_ref and determine if it's an Array or Hash */ |
| 1726 | 100 | T | F | if (!SvROK(data_ref)) |
| 1727 | 100 | T | F | croak("oneway_test: first argument must be a hash or array reference"); |
| 1730 | 100 | T | F | if (SvTYPE(rv) == SVt_PVHV) { |
| 1738 | 100 | T | F | /* MODE 3 – Array of Arrays (AoA) */ |
| 1742 | 100 | T | F | k = (size_t)av_len(in_av) + 1; |
| 1747 | 100 | T | F | /* first pass: sizes, total_n, and generate index names */ |
| 1750 | 50 | T | F | if (!val || !*val || !SvROK(*val) || SvTYPE(SvRV(*val)) != SVt_PVAV) |
| 1762 | 100 | T | F | memcpy(gnames[g], buf, klen + 1); |
| 1767 | 100 | T | F | for (size_t g = 0; g < k; g++) { |
| 1781 | 100 | T | F | factor_name = rhs; /* use the actual factor variable name */ |
| 1782 | 100 | T | F | SV **restrict resp_svp = hv_fetch(in_hv, lhs, (I32)strlen(lhs), 0); |
| 1805 | 50 | T | F | } else { |
| 50 | T | F | } else { |
| 1807 | 100 | T | F | k = (size_t)hv_iterinit(in_hv); |
| 50 | T | F | k = (size_t)hv_iterinit(in_hv); |
| 1812 | 50 | T | F | /* first pass: sizes, total_n, and group name strings */ |
| 1822 | 100 | T | F | croak("oneway_test: group '%s' has fewer than 2 observations", |
| 1823 | 50 | T | F | HePV(he, PL_na)); |
| 1826 | 50 | T | F | /* save a copy of the key string */ |
| 1827 | 100 | T | F | STRLEN klen; |
| 1828 | 50 | T | F | const char *kstr = HePV(he, klen); |
| 1829 | 100 | T | F | gnames[g] = (char *)safemalloc(klen + 1); |
| 1830 | 50 | T | F | memcpy(gnames[g], kstr, klen + 1); |
| 1834 | 50 | T | F | /* second pass: fill flat in the same iteration order */ |
| 50 | T | F | /* second pass: fill flat in the same iteration order */ |
| 1838 | 100 | T | F | hv_iterinit(in_hv); |
| 50 | T | F | hv_iterinit(in_hv); |
| 1842 | 50 | T | F | for (IV i = 0; i < len; i++) { |
| 50 | T | F | for (IV i = 0; i < len; i++) { |
| 1845 | 100 | T | F | } |
| 50 | T | F | } |
| 1846 | 100 | T | F | } |
| 50 | T | F | } |
| 1855 | 100 | T | F | for (size_t i = 0; i < sizes[g]; i++) sum += flat[offset + i]; |
| 1857 | 50 | T | F | offset += sizes[g]; |
| 1861 | 50 | T | F | res = c_oneway_test(flat, sizes, k, var_equal); |
| 50 | T | F | res = c_oneway_test(flat, sizes, k, var_equal); |
| 1865 | 100 | T | F | /* ── build return hash ref |
| 50 | T | F | /* ── build return hash ref |
| 1871 | 100 | T | F | */ |
| 1873 | 50 | T | F | /* Group (factor) sub-hash */ |
| 50 | T | F | /* Group (factor) sub-hash */ |
| 50 | T | F | /* Group (factor) sub-hash */ |
| 1874 | 0 | T | F | { |
| 1877 | 100 | T | F | hv_stores(g_hv, "Sum Sq", newSVnv(res.ss_between)); |
| 1880 | 100 | T | F | hv_stores(g_hv, "Pr(>F)", newSVnv(res.p_value)); |
| 1886 | 50 | T | F | HV *restrict r_hv = newHV(); |
| 1888 | 50 | T | F | hv_stores(r_hv, "Sum Sq", newSVnv(res.ss_within)); |
| 50 | T | F | hv_stores(r_hv, "Sum Sq", newSVnv(res.ss_within)); |
| 50 | T | F | hv_stores(r_hv, "Sum Sq", newSVnv(res.ss_within)); |
| 50 | T | F | hv_stores(r_hv, "Sum Sq", newSVnv(res.ss_within)); |
| 1892 | 100 | T | F | /* group_stats sub-hash */ |
| 1894 | 50 | T | F | HV *restrict gs_hv = newHV(); |
| 50 | T | F | HV *restrict gs_hv = newHV(); |
| 50 | T | F | HV *restrict gs_hv = newHV(); |
| 50 | T | F | HV *restrict gs_hv = newHV(); |
| 1902 | 50 | T | F | } |
| 1905 | 50 | T | F | hv_stores(ret_hv, "group_stats", newRV_noinc((SV *)gs_hv)); |
| 100 | T | F | hv_stores(ret_hv, "group_stats", newRV_noinc((SV *)gs_hv)); |
| 1909 | 100 | T | F | for (size_t g = 0; g < k; g++) Safefree(gnames[g]); |
| 1910 | 100 | T | F | Safefree(gnames); |
| 50 | T | F | Safefree(gnames); |
| 1912 | 100 | T | F | /* freed here, after factor_name is no longer needed */ |
| 1914 | 50 | T | F | OUTPUT: |
| 50 | T | F | OUTPUT: |
| 1920 | 100 | T | F | SV *restrict x_sv = NULL, *restrict y_sv = NULL; |
| 1924 | 100 | T | F | |
| 1930 | 100 | T | F | // Check if second argument is an array (2-sample) or a string representing a CDF (1-sample) |
| 1935 | 100 | T | F | } else if (SvPOK(ST(arg_idx))) { |
| 1943 | 50 | T | F | SV *restrict val = ST(arg_idx + 1); |
| 1944 | 100 | T | F | if (strEQ(key, "x")) x_sv = val; |
| 1946 | 50 | T | F | else if (strEQ(key, "exact")) { |
| 50 | T | F | else if (strEQ(key, "exact")) { |
| 1953 | 100 | T | F | |
| 1961 | 100 | T | F | |
| 1963 | 50 | T | F | croak("ks_test: alternative must be 'two.sided', 'less', or 'greater'"); |
| 1966 | 50 | T | F | AV *restrict x_av = (AV*)SvRV(x_sv); |
| 1968 | 100 | T | F | if (nx == 0) croak("Not enough 'x' observations"); |
| 1970 | 50 | T | F | // Extract 'x' array to C-array |
| 50 | T | F | // Extract 'x' array to C-array |
| 1971 | 50 | T | F | double *restrict x_data = (double *)safemalloc(nx * sizeof(double)); |
| 1972 | 100 | T | F | size_t valid_nx = 0; |
| 100 | T | F | size_t valid_nx = 0; |
| 1973 | 100 | T | F | for (size_t i = 0; i < nx; i++) { |
| 1977 | 50 | T | F | } |
| 1978 | 50 | T | F | } |
| 1990 | 100 | T | F | for (size_t i = 0; i < ny; i++) { |
| 1995 | 100 | T | F | } |
| 1998 | 100 | T | F | Safefree(x_data); Safefree(y_data); |
| 2001 | 100 | T | F | |
| 50 | T | F | |
| 2003 | 100 | T | F | calc_2sample_stats(x_data, valid_nx, y_data, valid_ny, &d, &d_plus, &d_minus); |
| 2005 | 50 | T | F | // Map alternative to the correct statistic |
| 50 | T | F | // Map alternative to the correct statistic |
| 2010 | 100 | T | F | // Determine if exact or asymptotic |
| 2015 | 100 | T | F | |
| 2018 | 50 | T | F | double *restrict comb = (double *)safemalloc(total_n * sizeof(double)); |
| 2019 | 100 | T | F | for(size_t i=0; i<valid_nx; i++) comb[i] = x_data[i]; |
| 50 | T | F | for(size_t i=0; i<valid_nx; i++) comb[i] = x_data[i]; |
| 2023 | 0 | T | F | bool has_ties = FALSE; |
| 2025 | 0 | T | F | if(comb[i] == comb[i-1]) { has_ties = TRUE; break; } |
| 0 | T | F | if(comb[i] == comb[i-1]) { has_ties = TRUE; break; } |
| 2027 | 0 | T | F | Safefree(comb); |
| 2037 | 100 | T | F | method_desc = "Two-sample Kolmogorov-Smirnov test (asymptotic)"; |
| 2038 | 100 | T | F | double z = statistic * sqrt((double)(valid_nx * valid_ny) / (valid_nx + valid_ny)); |
| 2040 | 50 | T | F | p_value = K2l(z, 0, 1e-9); |
| 50 | T | F | p_value = K2l(z, 0, 1e-9); |
| 2045 | 100 | T | F | Safefree(y_data); |
| 2047 | 100 | T | F | const char *restrict dist = SvPV_nolen(y_sv); |
| 2048 | 50 | T | F | if (strEQ(dist, "pnorm")) { |
| 2050 | 0 | T | F | double max_d = 0.0, max_d_plus = 0.0, max_d_minus = 0.0; |
| 0 | T | F | double max_d = 0.0, max_d_plus = 0.0, max_d_minus = 0.0; |
| 2053 | 0 | T | F | double cdf_obs_high = (double)(i + 1) / valid_nx; |
| 0 | T | F | double cdf_obs_high = (double)(i + 1) / valid_nx; |
| 2054 | 0 | T | F | double cdf_theor = approx_pnorm(x_data[i]); |
| 2057 | 0 | T | F | if (diff1 > max_d_plus) max_d_plus = diff1; |
| 2073 | 100 | T | F | warn("exact 1-sample 1-sided KS test not implemented; using asymptotic"); |
| 2075 | 50 | T | F | p_value = exp(-2.0 * z * z); |
| 50 | T | F | p_value = exp(-2.0 * z * z); |
| 2077 | 50 | T | F | } else { |
| 50 | T | F | } else { |
| 2080 | 100 | T | F | if (is_two_sided) p_value = K2l(z, 0, 1e-6); |
| 100 | T | F | if (is_two_sided) p_value = K2l(z, 0, 1e-6); |
| 2081 | 50 | T | F | else p_value = exp(-2.0 * z * z); |
| 2084 | 0 | T | F | Safefree(x_data); |
| 2096 | 100 | T | F | hv_stores(res, "p_value", newSVnv(p_value)); |
| 50 | T | F | hv_stores(res, "p_value", newSVnv(p_value)); |
| 2099 | 50 | T | F | RETVAL = newRV_noinc((SV*)res); |
| 2102 | 100 | T | F | RETVAL |
| 50 | T | F | RETVAL |
| 2104 | 100 | T | F | SV* wilcox_test(...) |
| 2106 | 50 | T | F | { |
| 50 | T | F | { |
| 2110 | 100 | T | F | short int exact = -1; |
| 2112 | 50 | T | F | int arg_idx = 0; |
| 50 | T | F | int arg_idx = 0; |
| 2116 | 100 | T | F | arg_idx++; |
| 2123 | 100 | T | F | // Ensure the remaining arguments form complete key-value pairs |
| 2128 | 100 | T | F | for (; arg_idx < items; arg_idx += 2) { |
| 2132 | 100 | T | F | else if (strEQ(key, "y")) y_sv = val; |
| 50 | T | F | else if (strEQ(key, "y")) y_sv = val; |
| 2135 | 0 | T | F | else if (strEQ(key, "mu")) mu = SvNV(val); |
| 2137 | 0 | T | F | if (!SvOK(val)) exact = -1; |
| 0 | T | F | if (!SvOK(val)) exact = -1; |
| 2139 | 0 | T | F | } |
| 2149 | 100 | T | F | |
| 2150 | 100 | T | F | AV *restrict y_av = NULL; |
| 2152 | 50 | T | F | if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) { |
| 50 | T | F | if (y_sv && SvROK(y_sv) && SvTYPE(SvRV(y_sv)) == SVt_PVAV) { |
| 2157 | 100 | T | F | const char *restrict method_desc = ""; |
| 2160 | 50 | T | F | if (ny > 0 && !paired) { |
| 50 | T | F | if (ny > 0 && !paired) { |
| 2161 | 100 | T | F | RankInfo *restrict ri = (RankInfo *)safemalloc((nx + ny) * sizeof(RankInfo)); |
| 2162 | 50 | T | F | size_t valid_nx = 0, valid_ny = 0; |
| 2163 | 0 | T | F | for (size_t i = 0; i < nx; i++) { |
| 2164 | 0 | T | F | SV**restrict el = av_fetch(x_av, i, 0); |
| 0 | T | F | SV**restrict el = av_fetch(x_av, i, 0); |
| 2165 | 0 | T | F | if (el && SvOK(*el) && looks_like_number(*el)) { |
| 2168 | 0 | T | F | valid_nx++; |
| 2182 | 100 | T | F | bool has_ties = 0; |
| 2184 | 50 | T | F | double w_rank_sum = 0.0; |
| 50 | T | F | double w_rank_sum = 0.0; |
| 2185 | 50 | T | F | for (size_t i = 0; i < total_n; i++) if (ri[i].idx == 1) w_rank_sum += ri[i].rank; |
| 2186 | 100 | T | F | statistic = w_rank_sum - (double)valid_nx * (valid_nx + 1.0) / 2.0; |
| 50 | T | F | statistic = w_rank_sum - (double)valid_nx * (valid_nx + 1.0) / 2.0; |
| 2187 | 50 | T | F | |
| 2190 | 0 | T | F | else use_exact = (valid_nx < 50 && valid_ny < 50 && !has_ties); |
| 2199 | 100 | T | F | double p_greater = 1.0 - exact_pwilcox(statistic - 1.0, valid_nx, valid_ny); |
| 50 | T | F | double p_greater = 1.0 - exact_pwilcox(statistic - 1.0, valid_nx, valid_ny); |
| 2203 | 50 | T | F | else { |
| 2204 | 100 | T | F | double p = (p_less < p_greater) ? p_less : p_greater; |
| 2221 | 50 | T | F | if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z); |
| 50 | T | F | if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z); |
| 50 | T | F | if (strcmp(alternative, "less") == 0) p_value = approx_pnorm(z); |
| 2226 | 50 | T | F | } else { // --- ONE SAMPLE / PAIRED --- |
| 2227 | 50 | T | F | if (paired && (!y_av || nx != ny)) croak("'x' and 'y' must have the same length for paired test"); |
| 2230 | 50 | T | F | bool has_zeroes = FALSE; |
| 2235 | 100 | T | F | |
| 2239 | 50 | T | F | double dy = SvNV(*y_el); |
| 100 | T | F | double dy = SvNV(*y_el); |
| 2241 | 50 | T | F | if (d == 0.0) has_zeroes = TRUE; // Drop exact zeroes |
| 100 | T | F | if (d == 0.0) has_zeroes = TRUE; // Drop exact zeroes |
| 2245 | 100 | T | F | if (d == 0.0) has_zeroes = TRUE; |
| 2248 | 50 | T | F | } |
| 2249 | 50 | T | F | if (n_nz == 0) { |
| 100 | T | F | if (n_nz == 0) { |
| 2251 | 50 | T | F | croak("not enough (non-missing) observations"); |
| 2254 | 50 | T | F | for (size_t i = 0; i < n_nz; i++) { |
| 50 | T | F | for (size_t i = 0; i < n_nz; i++) { |
| 50 | T | F | for (size_t i = 0; i < n_nz; i++) { |
| 2259 | 100 | T | F | double tie_adj = rank_and_count_ties(ri, n_nz, &has_ties); |
| 2261 | 100 | T | F | for (size_t i = 0; i < n_nz; i++) { |
| 2262 | 100 | T | F | if (ri[i].idx) statistic += ri[i].rank; |
| 2263 | 100 | T | F | } |
| 100 | T | F | } |
| 100 | T | F | } |
| 2266 | 100 | T | F | else use_exact = (n_nz < 50 && !has_ties); |
| 2269 | 100 | T | F | use_exact = FALSE; |
| 2272 | 100 | T | F | warn("cannot compute exact p-value with zeroes; falling back to approximation"); |
| 50 | T | F | warn("cannot compute exact p-value with zeroes; falling back to approximation"); |
| 50 | T | F | warn("cannot compute exact p-value with zeroes; falling back to approximation"); |
| 100 | T | F | warn("cannot compute exact p-value with zeroes; falling back to approximation"); |
| 2281 | 100 | T | F | else if (strcmp(alternative, "greater") == 0) p_value = p_greater; |
| 2290 | 50 | T | F | double z = statistic - exp; |
| 2294 | 50 | T | F | else if (strcmp(alternative, "greater") == 0) CORRECTION = 0.5; |
| 2295 | 50 | T | F | else if (strcmp(alternative, "less") == 0) CORRECTION = -0.5; |
| 2298 | 50 | T | F | |
| 2310 | 100 | T | F | hv_stores(res, "alternative", newSVpv(alternative, 0)); |
| 2312 | 50 | T | F | } |
| 2316 | 50 | T | F | SV* chisq_test(data_ref) |
| 2317 | 50 | T | F | SV* data_ref; |
| 2320 | 50 | T | F | // 1. Input Validation (mimics: die 'Input must be an array reference') |
| 2330 | 50 | T | F | is_2d = 1; |
| 2342 | 50 | T | F | double *restrict row_sum = (double*)safemalloc(r * sizeof(double)); |
| 50 | T | F | double *restrict row_sum = (double*)safemalloc(r * sizeof(double)); |
| 2345 | 50 | T | F | for(unsigned int j=0; j<c; j++) col_sum[j] = 0.0; |
| 50 | T | F | for(unsigned int j=0; j<c; j++) col_sum[j] = 0.0; |
| 2350 | 100 | T | F | SV**restrict val_sv = av_fetch(row, j, 0); |
| 2351 | 100 | T | F | double val = SvNV(*val_sv); |
| 2352 | 50 | T | F | row_sum[i] += val; |
| 2361 | 50 | T | F | for (unsigned int j = 0; j < c; j++) { |
| 2372 | 100 | T | F | } else { |
| 2377 | 50 | T | F | } |
| 50 | T | F | } |
| 50 | T | F | } |
| 2378 | 50 | T | F | safefree(row_sum); safefree(col_sum); |
| 50 | T | F | safefree(row_sum); safefree(col_sum); |
| 50 | T | F | safefree(row_sum); safefree(col_sum); |
| 2381 | 50 | T | F | for (unsigned int j = 0; j < c; j++) { |
| 50 | T | F | for (unsigned int j = 0; j < c; j++) { |
| 2389 | 50 | T | F | av_push(expected_av, newSVnv(E)); |
| 2395 | 100 | T | F | |
| 2397 | 100 | T | F | HV*restrict results = newHV(); |
| 2398 | 100 | T | F | |
| 2406 | 100 | T | F | hv_store(parameter_hv, "df", 2, newSViv(df), 0); |
| 2415 | 100 | T | F | // 'observed' => data_ref (Increment ref count since hv_store consumes ownership) |
| 2427 | 100 | T | F | } |
| 2483 | 50 | T | F | if (!data_sv || !SvROK(data_sv)) { |
| 2485 | 100 | T | F | } |
| 2488 | 100 | T | F | croak("write_table: 'data' must be a HASH or ARRAY reference\n"); |
| 2489 | 100 | T | F | } |
| 2490 | 50 | T | F | |
| 2493 | 50 | T | F | |
| 2494 | 50 | T | F | // NEW: Auto-detect separator from file extension if not overridden |
| 50 | T | F | // NEW: Auto-detect separator from file extension if not overridden |
| 2498 | 100 | T | F | const char *restrict ext = file + file_len - 4; |
| 50 | T | F | const char *restrict ext = file + file_len - 4; |
| 2506 | 100 | T | F | |
| 100 | T | F | |
| 50 | T | F | |
| 2510 | 50 | T | F | } |
| 2514 | 50 | T | F | AV *restrict rows_av = NULL; |
| 2515 | 50 | T | F | |
| 2518 | 100 | T | F | HV *restrict hv = (HV*)data_ref; |
| 2519 | 50 | T | F | if (hv_iterinit(hv) == 0) XSRETURN_EMPTY; |
| 2523 | 50 | T | F | if (!first_val || !SvROK(first_val)) { |
| 50 | T | F | if (!first_val || !SvROK(first_val)) { |
| 2528 | 50 | T | F | croak("write_table: Data values must be either all HASHes or all ARRAYs\n"); |
| 2531 | 0 | T | F | is_hoa = (first_type == SVt_PVAV); |
| 0 | T | F | is_hoa = (first_type == SVt_PVAV); |
| 2532 | 0 | T | F | hv_iterinit(hv); |
| 0 | T | F | hv_iterinit(hv); |
| 2541 | 50 | T | F | hv_iterinit(hv); |
| 0 | T | F | hv_iterinit(hv); |
| 2547 | 100 | T | F | AV *restrict av = (AV*)data_ref; |
| 2549 | 100 | T | F | SV **restrict first_ptr = av_fetch(av, 0, 0); |
| 2550 | 50 | T | F | if (!first_ptr || !*first_ptr || !SvROK(*first_ptr) || SvTYPE(SvRV(*first_ptr)) != SVt_PVHV) { |
| 2552 | 50 | T | F | } |
| 2558 | 50 | T | F | } |
| 2560 | 50 | T | F | is_aoh = 1; |
| 2562 | 50 | T | F | |
| 2564 | 50 | T | F | if (!fh) croak("write_table: Could not open '%s' for writing", file); |
| 100 | T | F | if (!fh) croak("write_table: Could not open '%s' for writing", file); |
| 2567 | 50 | T | F | bool inc_rownames = (row_names_sv && SvTRUE(row_names_sv)) ? 1 : 0; |
| 2568 | 100 | T | F | const char *restrict rownames_col = NULL; |
| 2572 | 50 | T | F | if (col_names_sv && SvOK(col_names_sv)) { |
| 50 | T | F | if (col_names_sv && SvOK(col_names_sv)) { |
| 2574 | 50 | T | F | for(size_t i=0; i<=av_len(c_av); i++) { |
| 50 | T | F | for(size_t i=0; i<=av_len(c_av); i++) { |
| 2576 | 100 | T | F | if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c)); |
| 2584 | 0 | T | F | hv_iterinit(inner); |
| 2587 | 0 | T | F | hv_store_ent(col_map, hv_iterkeysv(inner_entry), newSViv(1), 0); |
| 0 | T | F | hv_store_ent(col_map, hv_iterkeysv(inner_entry), newSViv(1), 0); |
| 2588 | 0 | T | F | } |
| 2590 | 0 | T | F | unsigned num_cols = hv_iterinit(col_map); |
| 0 | T | F | unsigned num_cols = hv_iterinit(col_map); |
| 0 | T | F | unsigned num_cols = hv_iterinit(col_map); |
| 2595 | 0 | T | F | } |
| 2603 | 100 | T | F | |
| 2604 | 50 | T | F | size_t h_idx = 0; |
| 2609 | 100 | T | F | } |
| 2612 | 100 | T | F | |
| 2614 | 50 | T | F | const char **restrict row_array = safemalloc(num_rows * sizeof(char*)); |
| 2615 | 100 | T | F | for(size_t i=0; i<num_rows; i++) { |
| 2617 | 50 | T | F | } |
| 2619 | 100 | T | F | |
| 2620 | 100 | T | F | HV *restrict data_hv = (HV*)data_ref; |
| 2622 | 100 | T | F | |
| 2623 | 50 | T | F | for(size_t i=0; i<num_rows; i++) { |
| 0 | T | F | for(size_t i=0; i<num_rows; i++) { |
| 2629 | 50 | T | F | |
| 2630 | 100 | T | F | for(size_t j=0; j<num_headers; j++) { |
| 2631 | 100 | T | F | SV**restrict h_ptr = av_fetch(headers_av, j, 0); |
| 2632 | 50 | T | F | const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : ""; |
| 2637 | 100 | T | F | safefree(row_array); safefree(row_data); |
| 2638 | 50 | T | F | if (headers_av) SvREFCNT_dec(headers_av); |
| 2649 | 100 | T | F | safefree(row_array); safefree(row_data); |
| 2660 | 50 | T | F | } |
| 50 | T | F | } |
| 2661 | 50 | T | F | |
| 2664 | 100 | T | F | for(size_t i=0; i<=av_len(c_av); i++) { |
| 2666 | 50 | T | F | if(c && SvOK(*c)) av_push(headers_av, newSVsv(*c)); |
| 2670 | 100 | T | F | const char **restrict col_array = safemalloc(num_cols * sizeof(char*)); |
| 2671 | 100 | T | F | for(unsigned int i=0; i<num_cols; i++) { |
| 2673 | 100 | T | F | col_array[i] = SvPV_nolen(hv_iterkeysv(ce)); |
| 2675 | 50 | T | F | qsort(col_array, num_cols, sizeof(char*), cmp_string_wt); |
| 2676 | 100 | T | F | for(unsigned i=0; i<num_cols; i++) av_push(headers_av, newSVpv(col_array[i], 0)); |
| 2681 | 50 | T | F | rownames_col = SvPV_nolen(row_names_sv); |
| 2684 | 50 | T | F | for(size_t i=0; i<=av_len(headers_av); i++) { |
| 2686 | 100 | T | F | if (!h_ptr || !*h_ptr) continue; |
| 2692 | 50 | T | F | SvREFCNT_dec(headers_av); |
| 2693 | 0 | T | F | headers_av = filtered_headers; |
| 2702 | 100 | T | F | } |
| 2705 | 100 | T | F | const char **restrict row_data = safemalloc((num_headers + 1) * sizeof(char*)); |
| 2707 | 100 | T | F | size_t d_idx = 0; |
| 2708 | 100 | T | F | if (inc_rownames) { |
| 2709 | 50 | T | F | if (rownames_col) { |
| 50 | T | F | if (rownames_col) { |
| 2713 | 100 | T | F | SV **restrict rn_val_ptr = av_fetch(rn_arr, i, 0); |
| 2714 | 50 | T | F | if (rn_val_ptr && SvOK(*rn_val_ptr)) { |
| 2723 | 50 | T | F | row_data[d_idx++] = undef_val; |
| 2724 | 100 | T | F | } |
| 2725 | 100 | T | F | } else { |
| 2728 | 50 | T | F | } else { |
| 2737 | 100 | T | F | SV **restrict arr_ptr = hv_fetch(data_hv, col_name, strlen(col_name), 0); |
| 100 | T | F | SV **restrict arr_ptr = hv_fetch(data_hv, col_name, strlen(col_name), 0); |
| 2738 | 100 | T | F | if (arr_ptr && SvROK(*arr_ptr)) { |
| 2740 | 100 | T | F | SV **restrict cell_ptr = av_fetch(arr, i, 0); |
| 2743 | 100 | T | F | PerlIO_close(fh); |
| 2747 | 100 | T | F | } |
| 2748 | 50 | T | F | row_data[d_idx++] = SvPV_nolen(*cell_ptr); |
| 2750 | 50 | T | F | row_data[d_idx++] = undef_val; |
| 100 | T | F | row_data[d_idx++] = undef_val; |
| 2756 | 100 | T | F | print_string_row(fh, row_data, d_idx, sep); |
| 2758 | 100 | T | F | } |
| 2760 | 50 | T | F | } else if (is_aoh) {// ----- Array of Hashes ----- |
| 100 | T | F | } else if (is_aoh) {// ----- Array of Hashes ----- |
| 2762 | 100 | T | F | size_t num_rows = av_len(data_av) + 1; |
| 2765 | 50 | T | F | for(size_t i=0; i<=av_len(c_av); i++) { |
| 2766 | 50 | T | F | SV **restrict c = av_fetch(c_av, i, 0); |
| 2769 | 100 | T | F | } else { |
| 2770 | 50 | T | F | HV *restrict col_map = newHV(); |
| 2780 | 100 | T | F | } |
| 100 | T | F | } |
| 50 | T | F | } |
| 2785 | 100 | T | F | HE *restrict ce = hv_iternext(col_map); |
| 2788 | 100 | T | F | qsort(col_array, num_cols, sizeof(char*), cmp_string_wt); |
| 2792 | 100 | T | F | } |
| 2795 | 100 | T | F | AV *restrict filtered_headers = newAV(); |
| 100 | T | F | AV *restrict filtered_headers = newAV(); |
| 2796 | 100 | T | F | for(size_t i=0; i<=av_len(headers_av); i++) { |
| 2797 | 100 | T | F | SV**restrict h_ptr = av_fetch(headers_av, i, 0); |
| 2798 | 50 | T | F | if (!h_ptr || !*h_ptr) continue; |
| 2799 | 100 | T | F | SV *restrict h_sv = *h_ptr; |
| 2801 | 100 | T | F | av_push(filtered_headers, newSVsv(h_sv)); |
| 2807 | 100 | T | F | size_t num_headers = av_len(headers_av) + 1; |
| 2808 | 100 | T | F | const char **restrict header_row = safemalloc((num_headers + 1) * sizeof(char*)); |
| 2809 | 100 | T | F | size_t h_idx = 0; |
| 2810 | 50 | T | F | if (inc_rownames) header_row[h_idx++] = ""; |
| 2818 | 100 | T | F | for(size_t i=0; i<num_rows; i++) { |
| 2821 | 50 | T | F | HV *restrict row_hv = (row_ptr && SvROK(*row_ptr)) ? (HV*)SvRV(*row_ptr) : NULL; |
| 2827 | 100 | T | F | PerlIO_close(fh); |
| 50 | T | F | PerlIO_close(fh); |
| 2828 | 100 | T | F | safefree(row_data); |
| 2830 | 100 | T | F | croak("write_table: Cannot write nested reference types to table\n"); |
| 2833 | 100 | T | F | } else { |
| 2834 | 50 | T | F | row_data[d_idx++] = undef_val; |
| 2836 | 100 | T | F | } else { |
| 2845 | 100 | T | F | const char *restrict col_name = (h_ptr && SvOK(*h_ptr)) ? SvPV_nolen(*h_ptr) : ""; |
| 2850 | 50 | T | F | safefree(row_data); |
| 2853 | 0 | T | F | } |
| 2854 | 0 | T | F | row_data[d_idx++] = SvPV_nolen(*cell_ptr); |
| 2858 | 100 | T | F | } |
| 2862 | 100 | T | F | safefree(row_data); |
| 2863 | 100 | T | F | } |
| 2885 | 100 | T | F | data = newAV(); |
| 2887 | 100 | T | F | sep_len = sep_str ? strlen(sep_str) : 0; |
| 2889 | 100 | T | F | |
| 2891 | 100 | T | F | if (!fp) { |
| 2898 | 100 | T | F | size_t len = SvCUR(line_sv); |
| 2908 | 50 | T | F | bool is_empty = 1; |
| 50 | T | F | bool is_empty = 1; |
| 2920 | 100 | T | F | for (size_t i = 0; i < len; i++) { |
| 2924 | 100 | T | F | if (in_quotes && (i + 1 < len) && line[i+1] == '"') { |
| 2925 | 100 | T | F | sv_catpvn(field, "\"", 1); |
| 2926 | 50 | T | F | i++; // Skip the escaped second quote |
| 2927 | 100 | T | F | } else if (in_quotes) { |
| 50 | T | F | } else if (in_quotes) { |
| 2928 | 50 | T | F | in_quotes = 0; // Close quotes |
| 2941 | 50 | T | F | } |
| 50 | T | F | } |
| 50 | T | F | } |
| 2942 | 50 | T | F | if (in_quotes) { |
| 50 | T | F | if (in_quotes) { |
| 50 | T | F | if (in_quotes) { |
| 2950 | 50 | T | F | // If a callback is provided, invoke it in a streaming fashion |
| 2956 | 100 | T | F | XPUSHs(sv_2mortal(newRV_inc((SV*)current_row))); |
| 2960 | 50 | T | F | LEAVE; |
| 100 | T | F | LEAVE; |
| 50 | T | F | LEAVE; |
| 2961 | 50 | T | F | SvREFCNT_dec(current_row); // Frees the row from C memory if Perl didn't keep it |
| 100 | T | F | SvREFCNT_dec(current_row); // Frees the row from C memory if Perl didn't keep it |
| 50 | T | F | SvREFCNT_dec(current_row); // Frees the row from C memory if Perl didn't keep it |
| 2964 | 100 | T | F | } |
| 100 | T | F | } |
| 2971 | 50 | T | F | if (in_quotes) { |
| 2977 | 100 | T | F | PUSHMARK(SP); |
| 2980 | 100 | T | F | call_sv(callback, G_DISCARD); |
| 2989 | 50 | T | F | SvREFCNT_dec(field); |
| 50 | T | F | SvREFCNT_dec(field); |
| 3003 | 100 | T | F | if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV) { |
| 3005 | 100 | T | F | } |
| 3006 | 100 | T | F | if (!SvROK(y_sv) || SvTYPE(SvRV(y_sv)) != SVt_PVAV) { |
| 3010 | 50 | T | F | // 2. Validate method argument |
| 0 | T | F | // 2. Validate method argument |
| 3011 | 50 | T | F | if (strcmp(method, "pearson") != 0 && |
| 3012 | 50 | T | F | strcmp(method, "spearman") != 0 && |
| 3013 | 100 | T | F | strcmp(method, "kendall") != 0) { |
| 3018 | 50 | T | F | AV *restrict y_av = (AV*)SvRV(y_sv); |
| 3020 | 50 | T | F | size_t ny = av_len(y_av) + 1; |
| 50 | T | F | size_t ny = av_len(y_av) + 1; |
| 3024 | 50 | T | F | (unsigned long)nx, (unsigned long)ny); |
| 0 | T | F | (unsigned long)nx, (unsigned long)ny); |
| 3025 | 50 | T | F | } |
| 50 | T | F | } |
| 3030 | 50 | T | F | double *restrict y_val = (double*)safemalloc(nx * sizeof(double)); |
| 50 | T | F | double *restrict y_val = (double*)safemalloc(nx * sizeof(double)); |
| 3032 | 50 | T | F | |
| 3040 | 0 | T | F | |
| 0 | T | F | |
| 3043 | 0 | T | F | x_val[n] = xv; |
| 3045 | 0 | T | F | n++; |
| 3051 | 50 | T | F | Safefree(x_val); Safefree(y_val); |
| 3059 | 100 | T | F | for (size_t j = 0; j < n; j++) { |
| 3068 | 50 | T | F | // Spearman: Rank the data first, then run standard covariance |
| 50 | T | F | // Spearman: Rank the data first, then run standard covariance |
| 3072 | 100 | T | F | rank_data(x_val, rx, n); |
| 3079 | 100 | T | F | cov_sum += dx * (ry[i] - mean_y); |
| 3080 | 50 | T | F | } |
| 50 | T | F | } |
| 3085 | 50 | T | F | double dx = x_val[i] - mean_x; |
| 0 | T | F | double dx = x_val[i] - mean_x; |
| 3086 | 50 | T | F | mean_x += dx / (i + 1); |
| 50 | T | F | mean_x += dx / (i + 1); |
| 3091 | 50 | T | F | } |
| 3096 | 0 | T | F | Safefree(x_val); Safefree(y_val); |
| 3113 | 100 | T | F | char **restrict terms = NULL, **restrict uniq_terms = NULL, **restrict exp_terms = NULL; |
| 3134 | 50 | T | F | HV *restrict res_hv, *restrict coef_hv, *restrict fitted_hv, *restrict resid_hv, *restrict summary_hv; |
| 50 | T | F | HV *restrict res_hv, *restrict coef_hv, *restrict fitted_hv, *restrict resid_hv, *restrict summary_hv; |
| 3141 | 50 | T | F | const char *restrict key = SvPV_nolen(ST(i_arg)); |
| 3144 | 100 | T | F | else if (strEQ(key, "data")) data_sv = val; |
| 3146 | 50 | T | F | else croak("glm: unknown argument '%s'", key); |
| 50 | T | F | else croak("glm: unknown argument '%s'", key); |
| 3148 | 50 | T | F | if (!formula) croak("glm: formula is required"); |
| 3156 | 50 | T | F | Newx(terms, term_cap, char*); Newx(uniq_terms, term_cap, char*); |
| 50 | T | F | Newx(terms, term_cap, char*); Newx(uniq_terms, term_cap, char*); |
| 3163 | 100 | T | F | |
| 3166 | 50 | T | F | *tilde = '\0'; |
| 3173 | 50 | T | F | while (chunk != NULL) { |
| 3177 | 0 | T | F | } |
| 3183 | 50 | T | F | if (star) { |
| 3184 | 50 | T | F | *star = '\0'; |
| 3185 | 100 | T | F | char *restrict left = chunk; char *restrict right = star + 1; |
| 3193 | 50 | T | F | snprintf(terms[num_terms++], inter_len, "%s:%s", left, right); |
| 100 | T | F | snprintf(terms[num_terms++], inter_len, "%s:%s", left, right); |
| 3195 | 100 | T | F | char *restrict c_chunk = strchr(chunk, '^'); |
| 3203 | 100 | T | F | bool found = FALSE; |
| 3207 | 100 | T | F | if (!found) uniq_terms[num_uniq++] = savepv(terms[i]); |
| 3217 | 100 | T | F | if (entry) { |
| 3223 | 50 | T | F | for(i = 0; i < n; i++) { |
| 3246 | 50 | T | F | row_hashes[i] = (HV*)SvRV(*val); |
| 3247 | 50 | T | F | char buf[32]; snprintf(buf, sizeof(buf), "%lu", i + 1); |
| 3257 | 50 | T | F | // --- Categorical Expansion --- |
| 3267 | 100 | T | F | if (is_column_categorical(data_hoa, row_hashes, n, uniq_terms[j])) { |
| 3269 | 100 | T | F | Newx(levels, levels_cap, char*); |
| 50 | T | F | Newx(levels, levels_cap, char*); |
| 3272 | 100 | T | F | if (str_val) { |
| 3274 | 50 | T | F | for (size_t l = 0; l < num_levels; l++) { |
| 100 | T | F | for (size_t l = 0; l < num_levels; l++) { |
| 3276 | 100 | T | F | } |
| 50 | T | F | } |
| 3285 | 100 | T | F | for (size_t l1 = 0; l1 < num_levels - 1; l1++) { |
| 3287 | 100 | T | F | if (strcmp(levels[l1], levels[l2]) > 0) { |
| 100 | T | F | if (strcmp(levels[l1], levels[l2]) > 0) { |
| 3296 | 100 | T | F | Renew(dummy_base, exp_cap, char*); Renew(dummy_level, exp_cap, char*); |
| 3297 | 50 | T | F | } |
| 3308 | 100 | T | F | } |
| 3310 | 100 | T | F | exp_terms[p_exp] = savepv(uniq_terms[j]); is_dummy[p_exp] = FALSE; p_exp++; |
| 50 | T | F | exp_terms[p_exp] = savepv(uniq_terms[j]); is_dummy[p_exp] = FALSE; p_exp++; |
| 3313 | 100 | T | F | p = p_exp; |
| 3315 | 50 | T | F | Newx(X, n * p, double); Newx(Y, n, double); |
| 100 | T | F | Newx(X, n * p, double); Newx(Y, n, double); |
| 3317 | 100 | T | F | |
| 100 | T | F | |
| 3326 | 100 | T | F | if (strcmp(exp_terms[j], "Intercept") == 0) { |
| 3328 | 100 | T | F | } else if (is_dummy[j]) { |
| 100 | T | F | } else if (is_dummy[j]) { |
| 3337 | 100 | T | F | } |
| 3338 | 50 | T | F | } |
| 3353 | 50 | T | F | W = (double*)safemalloc(valid_n * sizeof(double)); Z = (double*)safemalloc(valid_n * sizeof(double)); |
| 3357 | 100 | T | F | for (i = 0; i < p; i++) { beta[i] = 0.0; beta_old[i] = 0.0; } |
| 3359 | 100 | T | F | double sum_y = 0.0; |
| 100 | T | F | double sum_y = 0.0; |
| 3361 | 100 | T | F | double mean_y = sum_y / valid_n; |
| 3366 | 100 | T | F | eta[i] = log(mu[i] / (1.0 - mu[i])); |
| 3371 | 50 | T | F | deviance_old += dev; |
| 3380 | 100 | T | F | if (is_binomial) { |
| 3383 | 100 | T | F | if (varmu < 1e-10) varmu = 1e-10; |
| 3386 | 50 | T | F | } else { |
| 3394 | 50 | T | F | double w = W[k], z = Z[k]; |
| 3398 | 50 | T | F | for (size_t j = 0; j < p; j++) XtWX[i * p + j] += xw * X[k * p + j]; |
| 3400 | 50 | T | F | } |
| 3404 | 100 | T | F | double sum = 0.0; |
| 3406 | 50 | T | F | beta[i] = sum; |
| 3422 | 50 | T | F | |
| 3423 | 100 | T | F | double dev = 0.0; |
| 3431 | 100 | T | F | deviance_new += res * res; |
| 3435 | 100 | T | F | if (!is_binomial || deviance_new <= deviance_old + 1e-7 || !isfinite(deviance_new)) { |
| 3436 | 100 | T | F | continue; |
| 3437 | 50 | T | F | } |
| 3442 | 100 | T | F | // Convergence Check |
| 100 | T | F | // Convergence Check |
| 3443 | 100 | T | F | if (fabs(deviance_new - deviance_old) / (0.1 + fabs(deviance_new)) < epsilon) { |
| 100 | T | F | if (fabs(deviance_new - deviance_old) / (0.1 + fabs(deviance_new)) < epsilon) { |
| 3446 | 50 | T | F | deviance_old = deviance_new; |
| 3448 | 100 | T | F | } |
| 3463 | 100 | T | F | if (is_binomial) { |
| 100 | T | F | if (is_binomial) { |
| 3468 | 100 | T | F | double diff = Y[i] - wtdmu; |
| 3472 | 50 | T | F | // --- AIC Calculation --- |
| 3476 | 100 | T | F | } else if (is_binomial) { |
| 3478 | 50 | T | F | } |
| 50 | T | F | } |
| 3481 | 100 | T | F | df_res = valid_n - final_rank; |
| 3482 | 100 | T | F | dispersion = is_binomial ? 1.0 : ((df_res > 0) ? (deviance_new / df_res) : NAN); |
| 3485 | 50 | T | F | if (is_binomial) { |
| 3492 | 50 | T | F | } |
| 3495 | 50 | T | F | Safefree(valid_row_names[i]); |
| 3497 | 50 | T | F | Safefree(valid_row_names); |
| 3499 | 50 | T | F | summary_hv = newHV(); terms_av = newAV(); |
| 50 | T | F | summary_hv = newHV(); terms_av = newAV(); |
| 3505 | 50 | T | F | if (aliased[j]) { |
| 0 | T | F | if (aliased[j]) { |
| 3509 | 50 | T | F | hv_store(row_hv, is_binomial ? "Pr(>|z|)" : "Pr(>|t|)", 8, newSVpv("NaN", 0), 0); |
| 3513 | 50 | T | F | double p_val = is_binomial ? 2.0 * (1.0 - approx_pnorm(fabs(val_stat))) : get_t_pvalue(val_stat, df_res, "two.sided"); |
| 3514 | 50 | T | F | |
| 3515 | 50 | T | F | hv_store(row_hv, "Estimate", 8, newSVnv(beta[j]), 0); |
| 3516 | 50 | T | F | hv_store(row_hv, "Std. Error", 10, newSVnv(se), 0); |
| 3520 | 100 | T | F | hv_store(summary_hv, exp_terms[j], strlen(exp_terms[j]), newRV_noinc((SV*)row_hv), 0); |
| 3533 | 100 | T | F | hv_store(res_hv, "iter", 4, newSVuv(iter > max_iter ? max_iter : iter), 0); |
| 3535 | 100 | T | F | hv_store(res_hv, "rank", 4, newSVuv(final_rank), 0); |
| 3563 | 50 | T | F | if (items < 2 || items % 2 != 0) |
| 100 | T | F | if (items < 2 || items % 2 != 0) |
| 50 | T | F | if (items < 2 || items % 2 != 0) |
| 3569 | 50 | T | F | const char *restrict method = "pearson"; |
| 3572 | 100 | T | F | bool continuity = 0; |
| 3576 | 100 | T | F | const char *restrict key = SvPV_nolen(ST(i)); |
| 3577 | 50 | T | F | SV *restrict val = ST(i + 1); |
| 3580 | 50 | T | F | else if (strEQ(key, "method")) method = SvPV_nolen(val); |
| 50 | T | F | else if (strEQ(key, "method")) method = SvPV_nolen(val); |
| 50 | T | F | else if (strEQ(key, "method")) method = SvPV_nolen(val); |
| 3584 | 50 | T | F | else croak("cor_test: unknown argument '%s'", key); |
| 3588 | 50 | T | F | double *restrict x, *restrict y; |
| 3590 | 100 | T | F | |
| 3592 | 50 | T | F | bool is_kendall = (strcmp(method, "kendall") == 0); |
| 50 | T | F | bool is_kendall = (strcmp(method, "kendall") == 0); |
| 3596 | 50 | T | F | if (!SvOK(x_ref) || !SvROK(x_ref) || SvTYPE(SvRV(x_ref)) != SVt_PVAV || |
| 3607 | 50 | T | F | x = safemalloc(n_raw * sizeof(double)); |
| 50 | T | F | x = safemalloc(n_raw * sizeof(double)); |
| 50 | T | F | x = safemalloc(n_raw * sizeof(double)); |
| 3611 | 100 | T | F | for (size_t i = 0; i < n_raw; i++) { |
| 3613 | 50 | T | F | SV **restrict y_val = av_fetch(y_av, i, 0); |
| 50 | T | F | SV **restrict y_val = av_fetch(y_av, i, 0); |
| 3614 | 50 | T | F | |
| 50 | T | F | |
| 3621 | 0 | T | F | y[n] = yv; |
| 3627 | 100 | T | F | Safefree(x); |
| 3630 | 100 | T | F | } |
| 3632 | 100 | T | F | if (is_pearson) { |
| 3634 | 100 | T | F | double mean_x = 0.0, mean_y = 0.0, M2_x = 0.0, M2_y = 0.0, cov = 0.0; |
| 3648 | 50 | T | F | // Confidence interval using Fisher's Z transform |
| 3672 | 100 | T | F | double denom = sqrt((double)(c + d + tie_x) * (double)(c + d + tie_y)); |
| 3674 | 100 | T | F | |
| 50 | T | F | |
| 3677 | 100 | T | F | |
| 3679 | 50 | T | F | if (!exact_sv || !SvOK(exact_sv)) { |
| 100 | T | F | if (!exact_sv || !SvOK(exact_sv)) { |
| 3686 | 100 | T | F | |
| 3693 | 100 | T | F | double var_S = n * (n - 1) * (2.0 * n + 5.0) / 18.0; |
| 3694 | 100 | T | F | double S = c - d; |
| 3704 | 100 | T | F | } |
| 3706 | 100 | T | F | } else if (is_spearman) { |
| 50 | T | F | } else if (is_spearman) { |
| 3709 | 100 | T | F | compute_ranks(x, rank_x, n); |
| 3711 | 50 | T | F | |
| 100 | T | F | |
| 3718 | 100 | T | F | mean_y += dy / (i + 1); |
| 3725 | 50 | T | F | // S = sum of squared rank differences (R's reported statistic) |
| 3726 | 100 | T | F | double S_stat = 0.0; |
| 3737 | 100 | T | F | break; |
| 3739 | 100 | T | F | } |
| 50 | T | F | } |
| 3742 | 100 | T | F | } else { |
| 3744 | 50 | T | F | } |
| 100 | T | F | } |
| 3754 | 100 | T | F | p_value = get_t_pvalue(statistic, (double)(n - 2), alternative); |
| 3764 | 100 | T | F | hv_stores(rhv, "p.value", newSVnv(p_value)); |
| 3765 | 50 | T | F | hv_stores(rhv, "statistic", newSVnv(statistic)); |
| 3777 | 100 | T | F | } |
| 3779 | 100 | T | F | RETVAL |
| 50 | T | F | RETVAL |
| 3782 | 100 | T | F | SV *data |
| 3784 | 50 | T | F | AV *restrict av; |
| 100 | T | F | AV *restrict av; |
| 3794 | 100 | T | F | n_raw = av_len(av) + 1; |
| 3804 | 100 | T | F | x[n] = val; |
| 3805 | 100 | T | F | mean += val; |
| 3821 | 50 | T | F | if (ssq == 0.0) { |
| 100 | T | F | if (ssq == 0.0) { |
| 50 | T | F | if (ssq == 0.0) { |
| 3827 | 50 | T | F | // --- Core AS R94 Algorithm: Weights and Statistic W --- |
| 100 | T | F | // --- Core AS R94 Algorithm: Weights and Statistic W --- |
| 50 | T | F | // --- Core AS R94 Algorithm: Weights and Statistic W --- |
| 3833 | 50 | T | F | // Exact P-value for n=3 |
| 3838 | 100 | T | F | Newx(m, n, double); |
| 3842 | 100 | T | F | sum_m2 += m[i] * m[i]; |
| 3843 | 100 | T | F | } |
| 3844 | 100 | T | F | double u = 1.0 / sqrt((double)n); |
| 3845 | 100 | T | F | double a_n = -2.706056*pow(u,5) + 4.434685*pow(u,4) - 2.071190*pow(u,3) - 0.147981*pow(u,2) + 0.221157*u + m[n-1]/sqrt(sum_m2); |
| 3846 | 100 | T | F | a[n-1] = a_n; |
| 3847 | 100 | T | F | a[0] = -a_n; |
| 3848 | 50 | T | F | if (n == 4 || n == 5) { |
| 3853 | 100 | T | F | } else { |
| 50 | T | F | } else { |
| 50 | T | F | } else { |
| 3857 | 50 | T | F | double eps = (sum_m2 - 2.0 * m[n-1]*m[n-1] - 2.0 * m[n-2]*m[n-2]) / (1.0 - 2.0 * a_n*a_n - 2.0 * a_n1*a_n1); |
| 3859 | 100 | T | F | a[i] = m[i] / sqrt(eps); |
| 50 | T | F | a[i] = m[i] / sqrt(eps); |
| 50 | T | F | a[i] = m[i] / sqrt(eps); |
| 3862 | 50 | T | F | for (size_t i = 0; i < n; i++) { |
| 100 | T | F | for (size_t i = 0; i < n; i++) { |
| 3867 | 100 | T | F | /* NOTE: p_val is declared in PREINIT above; |
| 3869 | 50 | T | F | */ |
| 50 | T | F | */ |
| 3875 | 100 | T | F | // if y reaches gamma the p-value is essentially zero |
| 50 | T | F | // if y reaches gamma the p-value is essentially zero |
| 3877 | 100 | T | F | double gamma = 0.459 * nn - 2.273; |
| 100 | T | F | double gamma = 0.459 * nn - 2.273; |
| 3878 | 100 | T | F | if (y >= gamma) { |
| 3880 | 100 | T | F | } else { |
| 100 | T | F | } else { |
| 3882 | 100 | T | F | double mu = 0.544 + nn * (-0.39978 + nn * ( 0.025054 - nn * 0.0006714)); |
| 3884 | 50 | T | F | double sigma = exp(sig_val); |
| 50 | T | F | double sigma = exp(sig_val); |
| 3890 | 100 | T | F | } else { |
| 3892 | 100 | T | F | double ln_n = log((double)n); |
| 3895 | 50 | T | F | double sig_val= -0.4803 + ln_n * (-0.082676 + ln_n * 0.0030302); |
| 50 | T | F | double sig_val= -0.4803 + ln_n * (-0.082676 + ln_n * 0.0030302); |
| 3896 | 50 | T | F | double sigma = exp(sig_val); |
| 50 | T | F | double sigma = exp(sig_val); |
| 3903 | 50 | T | F | |
| 3909 | 100 | T | F | hv_stores(ret_hash, "W", newSVnv(w)); |
| 3910 | 50 | T | F | hv_stores(ret_hash, "p_value", newSVnv(p_val)); |
| 0 | T | F | hv_stores(ret_hash, "p_value", newSVnv(p_val)); |
| 3919 | 50 | T | F | size_t count = 0; |
| 0 | T | F | size_t count = 0; |
| 3939 | 100 | T | F | } |
| 3943 | 100 | T | F | min_val = val; |
| 3966 | 100 | T | F | AV* restrict av = (AV*)SvRV(arg); |
| 50 | T | F | AV* restrict av = (AV*)SvRV(arg); |
| 3972 | 100 | T | F | if (first || val > max_val) { |
| 3978 | 100 | T | F | croak("max: undefined value at array ref index %zu (argument %zu)", j, i); |
| 3980 | 100 | T | F | } |
| 100 | T | F | } |
| 3981 | 100 | T | F | } else if (SvOK(arg)) { |
| 50 | T | F | } else if (SvOK(arg)) { |
| 3982 | 50 | T | F | double val = SvNV(arg); |
| 3986 | 50 | T | F | } |
| 3987 | 50 | T | F | count++; |
| 3989 | 100 | T | F | croak("max: undefined value at argument index %zu", i); |
| 3991 | 50 | T | F | } |
| 50 | T | F | } |
| 3997 | 100 | T | F | SV* runif(...) |
| 3998 | 100 | T | F | CODE: |
| 4000 | 100 | T | F | size_t n = 0; |
| 4002 | 100 | T | F | |
| 4004 | 100 | T | F | bool n_set = 0, min_set = 0, max_set = 0; |
| 4006 | 100 | T | F | unsigned int i = 0; |
| 4007 | 100 | T | F | |
| 4009 | 100 | T | F | croak("Usage: runif(n, [min=0], [max=1]) or runif(n => $n, ...)"); |
| 4011 | 100 | T | F | |
| 4013 | 100 | T | F | // 1. Check if the current argument is a string key for a named parameter |
| 4014 | 50 | T | F | if (i + 1 < items && SvPOK(ST(i))) { |
| 4016 | 100 | T | F | if (strEQ(key, "n")) { |
| 4018 | 100 | T | F | n_set = 1; |
| 4020 | 100 | T | F | continue; |
| 4021 | 50 | T | F | } else if (strEQ(key, "min")) { |
| 4023 | 100 | T | F | min_set = 1; |
| 4025 | 100 | T | F | continue; |
| 4027 | 100 | T | F | max = SvNV(ST(i+1)); |
| 4029 | 100 | T | F | i += 2; |
| 4030 | 100 | T | F | continue; |
| 4032 | 100 | T | F | } |
| 4034 | 50 | T | F | // 2. Fallback to positional parsing if it's not a recognized key |
| 4035 | 50 | T | F | if (!n_set) { |
| 4038 | 100 | T | F | } else if (!min_set) { |
| 4040 | 50 | T | F | min_set = 1; |
| 4045 | 100 | T | F | croak("Too many arguments or unrecognized parameter passed to runif()"); |
| 4049 | 100 | T | F | if (!n_set) { |
| 4054 | 100 | T | F | AV *restrict results = newAV(); |
| 4056 | 100 | T | F | av_extend(results, n - 1); |
| 4061 | 100 | T | F | if (max < min) { |
| 4063 | 100 | T | F | } else { |
| 4066 | 100 | T | F | av_push(results, newSVnv(r)); |
| 4070 | 100 | T | F | OUTPUT: |
| 4071 | 100 | T | F | RETVAL |
| 4077 | 100 | T | F | AUTO_SEED_PRNG(); |
| 4078 | 100 | T | F | if (items % 2 != 0) |
| 4079 | 50 | T | F | croak("Usage: rbinom(n => 10, size => 100, prob => 0.5)"); |
| 4083 | 50 | T | F | |
| 4084 | 0 | T | F | bool size_set = FALSE, prob_set = FALSE; |
| 4092 | 50 | T | F | else if (strEQ(key, "prob")) { prob = SvNV(val); prob_set = TRUE; } |
| 4093 | 100 | T | F | else croak("rbinom: unknown argument '%s'", key); |
| 4107 | 100 | T | F | |
| 4109 | 100 | T | F | } |
| 50 | T | F | } |
| 4112 | 100 | T | F | |
| 4114 | 50 | T | F | hist(SV* x_sv, ...) |
| 100 | T | F | hist(SV* x_sv, ...) |
| 4120 | 100 | T | F | |
| 4126 | 100 | T | F | double *restrict x; |
| 4129 | 50 | T | F | double min_val = DBL_MAX, max_val = -DBL_MAX; |
| 4132 | 100 | T | F | SV**restrict tv = av_fetch(x_av, i, 0); |
| 4134 | 100 | T | F | double val = SvNV(*tv); |
| 50 | T | F | double val = SvNV(*tv); |
| 4137 | 100 | T | F | if (val > max_val) max_val = val; |
| 4139 | 50 | T | F | } |
| 50 | T | F | } |
| 4146 | 50 | T | F | |
| 4156 | 100 | T | F | break; |
| 4163 | 100 | T | F | } |
| 4170 | 100 | T | F | Newx(density, n_bins, double); |
| 4171 | 100 | T | F | Newx(counts, n_bins, size_t); |
| 4172 | 100 | T | F | |
| 4177 | 50 | T | F | } |
| 50 | T | F | } |
| 4182 | 50 | T | F | // 6. Build Return HashRef |
| 4188 | 50 | T | F | for (size_t i = 0; i <= n_bins; i++) { |
| 100 | T | F | for (size_t i = 0; i <= n_bins; i++) { |
| 50 | T | F | for (size_t i = 0; i <= n_bins; i++) { |
| 4193 | 50 | T | F | av_push(av_density, newSVnv(density[i])); |
| 50 | T | F | av_push(av_density, newSVnv(density[i])); |
| 4194 | 50 | T | F | } |
| 4196 | 50 | T | F | hv_stores(res_hv, "breaks", newRV_noinc((SV*)av_breaks)); |
| 4197 | 50 | T | F | hv_stores(res_hv, "counts", newRV_noinc((SV*)av_counts)); |
| 4200 | 50 | T | F | |
| 50 | T | F | |
| 4202 | 50 | T | F | Safefree(x); Safefree(breaks); Safefree(mids); |
| 100 | T | F | Safefree(x); Safefree(breaks); Safefree(mids); |
| 50 | T | F | Safefree(x); Safefree(breaks); Safefree(mids); |
| 4208 | 100 | T | F | RETVAL |
| 50 | T | F | RETVAL |
| 4209 | 50 | T | F | |
| 4213 | 100 | T | F | SV *restrict x_sv = NULL; |
| 4217 | 50 | T | F | /* --- 1. Consume first positional arg as 'x' if it's an array ref --- */ |
| 4221 | 50 | T | F | } |
| 4222 | 50 | T | F | |
| 4224 | 100 | T | F | if ((items - arg_idx) % 2 != 0) |
| 4226 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 4228 | 100 | T | F | const char *restrict key = SvPV_nolen(ST(arg_idx)); |
| 4230 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 4239 | 50 | T | F | if (n_raw == 0) croak("quantile: 'x' is empty"); |
| 4244 | 50 | T | F | size_t n = 0; |
| 50 | T | F | size_t n = 0; |
| 50 | T | F | size_t n = 0; |
| 4248 | 50 | T | F | x[n++] = SvNV(*tv); |
| 4253 | 100 | T | F | croak("quantile: 'x' contains no valid numbers"); |
| 4255 | 50 | T | F | // --- Sort Data for Quantile Math --- |
| 50 | T | F | // --- Sort Data for Quantile Math --- |
| 50 | T | F | // --- Sort Data for Quantile Math --- |
| 4259 | 50 | T | F | unsigned int n_probs = 5; |
| 50 | T | F | unsigned int n_probs = 5; |
| 4260 | 50 | T | F | double *restrict probs; |
| 4261 | 100 | T | F | |
| 4263 | 50 | T | F | AV *restrict p_av = (AV*)SvRV(probs_sv); |
| 50 | T | F | AV *restrict p_av = (AV*)SvRV(probs_sv); |
| 50 | T | F | AV *restrict p_av = (AV*)SvRV(probs_sv); |
| 4270 | 50 | T | F | Safefree(x); Safefree(probs); |
| 4272 | 100 | T | F | } |
| 4273 | 50 | T | F | } |
| 4274 | 100 | T | F | } else { |
| 4278 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 4288 | 50 | T | F | q = x[n - 1]; /* Prevent out-of-bounds mapping */ |
| 50 | T | F | q = x[n - 1]; /* Prevent out-of-bounds mapping */ |
| 4292 | 50 | T | F | /* Continuous sample quantile interpolation (Type 7) */ |
| 4294 | 50 | T | F | unsigned int j = (unsigned int)h; /* floor via cast */ |
| 4295 | 100 | T | F | double gamma = h - j; |
| 4296 | 50 | T | F | q = (1.0 - gamma) * x[j] + gamma * x[j + 1]; |
| 4297 | 100 | T | F | } |
| 4301 | 50 | T | F | double pct = p * 100.0; |
| 50 | T | F | double pct = p * 100.0; |
| 50 | T | F | double pct = p * 100.0; |
| 4309 | 50 | T | F | hv_store(res_hv, key, strlen(key), newSVnv(q), 0); |
| 4316 | 50 | T | F | } |
| 4317 | 100 | T | F | OUTPUT: |
| 4321 | 50 | T | F | double mean(...) |
| 4324 | 0 | T | F | double total = 0; |
| 4325 | 0 | T | F | size_t count = 0; |
| 4326 | 0 | T | F | CODE: |
| 4328 | 0 | T | F | SV* restrict arg = ST(i); |
| 4330 | 0 | T | F | AV* restrict av = (AV*)SvRV(arg); |
| 4337 | 0 | T | F | } else { |
| 4338 | 0 | T | F | croak("mean: undefined value at array ref index %zu (argument %zu)", j, i); |
| 4341 | 0 | T | F | } else if (SvOK(arg)) { |
| 4345 | 100 | T | F | croak("mean: undefined value at argument index %zu", i); |
| 4346 | 100 | T | F | } |
| 4350 | 100 | T | F | OUTPUT: |
| 4354 | 100 | T | F | PROTOTYPE: @ |
| 4356 | 50 | T | F | HV *restrict counts; |
| 4357 | 100 | T | F | HV *restrict originals; |
| 4373 | 50 | T | F | if (tv && SvOK(*tv)) { |
| 4375 | 100 | T | F | const char *restrict key = SvPV(*tv, klen); |
| 100 | T | F | const char *restrict key = SvPV(*tv, klen); |
| 4380 | 50 | T | F | if (cnt > max_count) max_count = cnt; |
| 4382 | 50 | T | F | hv_store(originals, key, klen, newSVsv(*tv), 0); |
| 4387 | 50 | T | F | } |
| 50 | T | F | } |
| 100 | T | F | } |
| 4389 | 50 | T | F | STRLEN klen; |
| 50 | T | F | STRLEN klen; |
| 50 | T | F | STRLEN klen; |
| 0 | T | F | STRLEN klen; |
| 4391 | 0 | T | F | SV **restrict slot = hv_fetch(counts, key, klen, 1); |
| 4393 | 0 | T | F | size_t cnt = SvOK(*slot) ? SvIV(*slot) + 1 : 1; |
| 4402 | 100 | T | F | } |
| 4404 | 50 | T | F | if (arg_count == 0) |
| 4408 | 50 | T | F | while ((he = hv_iternext(counts))) { |
| 50 | T | F | while ((he = hv_iternext(counts))) { |
| 50 | T | F | while ((he = hv_iternext(counts))) { |
| 4410 | 50 | T | F | STRLEN klen; |
| 50 | T | F | STRLEN klen; |
| 50 | T | F | STRLEN klen; |
| 0 | T | F | STRLEN klen; |
| 4412 | 0 | T | F | SV **restrict orig = hv_fetch(originals, key, klen, 0); |
| 4414 | 0 | T | F | } |
| 4415 | 0 | T | F | } |
| 4426 | 100 | T | F | AV* restrict av = (AV*)SvRV(arg); |
| 4428 | 100 | T | F | for (size_t j = 0; j < len; j++) { |
| 50 | T | F | for (size_t j = 0; j < len; j++) { |
| 4430 | 50 | T | F | if (tv && SvOK(*tv)) { |
| 4432 | 50 | T | F | count++; |
| 50 | T | F | count++; |
| 50 | T | F | count++; |
| 4438 | 100 | T | F | total += SvNV(arg); |
| 4448 | 50 | T | F | |
| 50 | T | F | |
| 4454 | 100 | T | F | CODE: |
| 4460 | 100 | T | F | size_t len = av_len(av) + 1; |
| 4463 | 50 | T | F | if (tv && SvOK(*tv)) { |
| 4465 | 100 | T | F | double val = SvNV(*tv); |
| 4467 | 50 | T | F | mean += delta / count; |
| 50 | T | F | mean += delta / count; |
| 4470 | 50 | T | F | croak("sd: undefined value at array ref index %zu (argument %zu)", j, i); |
| 50 | T | F | croak("sd: undefined value at array ref index %zu (argument %zu)", j, i); |
| 4477 | 50 | T | F | mean += delta / count; |
| 4480 | 50 | T | F | croak("sd: undefined value at argument index %zu", i); |
| 4481 | 50 | T | F | } |
| 4487 | 100 | T | F | |
| 4494 | 100 | T | F | CODE: |
| 4496 | 50 | T | F | for (size_t i = 0; i < items; i++) { |
| 4503 | 50 | T | F | if (tv && SvOK(*tv)) { |
| 4512 | 100 | T | F | } |
| 4514 | 50 | T | F | count++; |
| 0 | T | F | count++; |
| 4517 | 0 | T | F | mean += delta / count; |
| 4519 | 0 | T | F | } else { |
| 0 | T | F | } else { |
| 4521 | 50 | T | F | } |
| 4525 | 50 | T | F | OUTPUT: |
| 4526 | 50 | T | F | RETVAL |
| 4527 | 100 | T | F | |
| 4529 | 50 | T | F | CODE: |
| 0 | T | F | CODE: |
| 4532 | 0 | T | F | SV*restrict y_sv = NULL; |
| 4534 | 0 | T | F | bool paired = FALSE, var_equal = FALSE; |
| 0 | T | F | bool paired = FALSE, var_equal = FALSE; |
| 4539 | 50 | T | F | // 1. Shift first positional argument as 'x' if it's an array reference |
| 4544 | 100 | T | F | |
| 4545 | 100 | T | F | // 2. Shift second positional argument as 'y' if it's an array reference |
| 4546 | 100 | T | F | if (arg_idx < items && SvROK(ST(arg_idx)) && SvTYPE(SvRV(ST(arg_idx))) == SVt_PVAV) { |
| 4551 | 100 | T | F | // Ensure the remaining arguments form complete key-value pairs |
| 4557 | 50 | T | F | for (; arg_idx < items; arg_idx += 2) { |
| 4558 | 100 | T | F | const char*restrict key = SvPV_nolen(ST(arg_idx)); |
| 4560 | 50 | T | F | |
| 4570 | 50 | T | F | |
| 4577 | 100 | T | F | AV*restrict y_av = NULL; |
| 4580 | 100 | T | F | |
| 4582 | 100 | T | F | croak("t_test: 'conf_level' must be between 0 and 1"); |
| 4585 | 50 | T | F | HV*restrict results = newHV(); |
| 4588 | 0 | T | F | double val = (tv && SvOK(*tv)) ? SvNV(*tv) : 0; |
| 4595 | 50 | T | F | |
| 100 | T | F | |
| 50 | T | F | |
| 4599 | 50 | T | F | if (paired && ny != nx) croak("t_test: Paired arrays must be same length"); |
| 4600 | 100 | T | F | double mean_y = 0.0, M2_y = 0.0, var_y; |
| 4604 | 50 | T | F | double delta = val - mean_y; |
| 0 | T | F | double delta = val - mean_y; |
| 4607 | 50 | T | F | } |
| 100 | T | F | } |
| 4609 | 50 | T | F | if (paired) { |
| 0 | T | F | if (paired) { |
| 4613 | 100 | T | F | SV**restrict dy_ptr = av_fetch(y_av, i, 0); |
| 50 | T | F | SV**restrict dy_ptr = av_fetch(y_av, i, 0); |
| 4621 | 100 | T | F | double var_d = M2_d / (nx - 1); |
| 4628 | 100 | T | F | } else if (var_equal) { |
| 4631 | 50 | T | F | cint_est = mean_x - mean_y; |
| 4632 | 50 | T | F | std_err = sqrt(pooled_var * (1.0 / nx + 1.0 / ny)); |
| 4675 | 50 | T | F | hv_store(results, "df", 2, newSVnv(df), 0); |
| 4677 | 100 | T | F | hv_store(results, "conf_int", 8, newRV_noinc((SV*)conf_int), 0); |
| 4680 | 100 | T | F | OUTPUT: |
| 4681 | 50 | T | F | RETVAL |
| 4684 | 100 | T | F | INIT: |
| 4685 | 100 | T | F | if (!SvROK(p_sv) || SvTYPE(SvRV(p_sv)) != SVt_PVAV) { |
| 100 | T | F | if (!SvROK(p_sv) || SvTYPE(SvRV(p_sv)) != SVt_PVAV) { |
| 4691 | 50 | T | F | if (n == 0) { |
| 4693 | 50 | T | F | } |
| 4695 | 50 | T | F | char meth[64]; |
| 4697 | 50 | T | F | for(unsigned short int i = 0; meth[i]; i++) meth[i] = tolower(meth[i]); |
| 100 | T | F | for(unsigned short int i = 0; meth[i]; i++) meth[i] = tolower(meth[i]); |
| 4700 | 50 | T | F | if (strstr(meth, "benjamini") && strstr(meth, "yekutieli")) strcpy(meth, "by"); |
| 4701 | 100 | T | F | if (strcmp(meth, "fdr") == 0) strcpy(meth, "bh"); |
| 4706 | 50 | T | F | Newx(adj, n, double); |
| 50 | T | F | Newx(adj, n, double); |
| 4708 | 50 | T | F | for (size_t i = 0; i < n; i++) { |
| 50 | T | F | for (size_t i = 0; i < n; i++) { |
| 4710 | 100 | T | F | arr[i].p = (tv && SvOK(*tv)) ? SvNV(*tv) : 1.0; |
| 4718 | 0 | T | F | double v = arr[i].p * n; |
| 4720 | 0 | T | F | } |
| 4721 | 0 | T | F | } else if (strcmp(meth, "holm") == 0) { |
| 4722 | 0 | T | F | double cummax = 0.0; |
| 4724 | 0 | T | F | double v = arr[i].p * (n - i); |
| 0 | T | F | double v = arr[i].p * (n - i); |
| 0 | T | F | double v = arr[i].p * (n - i); |
| 4729 | 0 | T | F | double cummin = 1.0; |
| 4740 | 100 | T | F | adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0; |
| 100 | T | F | adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0; |
| 50 | T | F | adj[arr[i].orig_idx] = (cummin < 1.0) ? cummin : 1.0; |
| 4744 | 100 | T | F | for (size_t i = 1; i <= n; i++) q += 1.0 / i; |
| 4745 | 100 | T | F | double cummin = 1.0; |
| 4746 | 50 | T | F | for (ssize_t i = n - 1; i >= 0; i--) { |
| 4758 | 100 | T | F | double temp = (n * arr[i].p) / (i + 1.0); |
| 4760 | 50 | T | F | min_val = temp; |
| 0 | T | F | min_val = temp; |
| 4762 | 0 | T | F | } |
| 0 | T | F | } |
| 0 | T | F | } |
| 0 | T | F | } |
| 4766 | 100 | T | F | q_arr[i] = min_val; |
| 50 | T | F | q_arr[i] = min_val; |
| 4767 | 50 | T | F | } |
| 0 | T | F | } |
| 0 | T | F | } |
| 4773 | 100 | T | F | for (size_t k = 1; k < i2_len; k++) { |
| 50 | T | F | for (size_t k = 1; k < i2_len; k++) { |
| 4774 | 0 | T | F | double temp_q1 = (j * arr[n_mj + 1 + k].p) / (2.0 + k); |
| 0 | T | F | double temp_q1 = (j * arr[n_mj + 1 + k].p) / (2.0 + k); |
| 0 | T | F | double temp_q1 = (j * arr[n_mj + 1 + k].p) / (2.0 + k); |
| 4780 | 100 | T | F | for (size_t i = 0; i <= n_mj; i++) { |
| 50 | T | F | for (size_t i = 0; i <= n_mj; i++) { |
| 0 | T | F | for (size_t i = 0; i <= n_mj; i++) { |
| 4786 | 100 | T | F | q_arr[n_mj + 1 + i] = q_arr[n_mj]; |
| 50 | T | F | q_arr[n_mj + 1 + i] = q_arr[n_mj]; |
| 0 | T | F | q_arr[n_mj + 1 + i] = q_arr[n_mj]; |
| 4790 | 100 | T | F | if (pa[i] < q_arr[i]) { |
| 50 | T | F | if (pa[i] < q_arr[i]) { |
| 4791 | 0 | T | F | pa[i] = q_arr[i]; |
| 0 | T | F | pa[i] = q_arr[i]; |
| 0 | T | F | pa[i] = q_arr[i]; |
| 4796 | 100 | T | F | for (size_t i = 0; i < n; i++) { |
| 4797 | 50 | T | F | double v = (pa[i] > arr[i].p) ? pa[i] : arr[i].p; |
| 0 | T | F | double v = (pa[i] > arr[i].p) ? pa[i] : arr[i].p; |
| 4798 | 50 | T | F | if (v > 1.0) v = 1.0; |
| 0 | T | F | if (v > 1.0) v = 1.0; |
| 4807 | 50 | T | F | Safefree(arr); Safefree(adj); |
| 4809 | 50 | T | F | } |
| 4811 | 50 | T | F | EXTEND(SP, n); |
| 50 | T | F | EXTEND(SP, n); |
| 4818 | 100 | T | F | double median(...) |
| 4819 | 100 | T | F | PROTOTYPE: @ |
| 4821 | 100 | T | F | size_t total_count = 0, k = 0; |
| 4823 | 50 | T | F | double median_val = 0.0; |
| 50 | T | F | double median_val = 0.0; |
| 4825 | 100 | T | F | /* Pass 1: Count valid elements — die immediately on any undef */ |
| 4827 | 50 | T | F | SV* restrict arg = ST(i); |
| 4828 | 100 | T | F | if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) { |
| 4838 | 50 | T | F | } |
| 4839 | 100 | T | F | } else if (SvOK(arg)) { |
| 4851 | 100 | T | F | for (size_t i = 0; i < items; i++) { |
| 4853 | 50 | T | F | if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) { |
| 4855 | 100 | T | F | size_t len = av_len(av) + 1; |
| 4856 | 50 | T | F | for (size_t j = 0; j < len; j++) { |
| 4861 | 100 | T | F | Safefree(nums); |
| 4866 | 50 | T | F | nums[k++] = SvNV(arg); |
| 0 | T | F | nums[k++] = SvNV(arg); |
| 4868 | 50 | T | F | Safefree(nums); |
| 50 | T | F | Safefree(nums); |
| 4876 | 50 | T | F | median_val = (nums[total_count / 2 - 1] + nums[total_count / 2]) / 2.0; |
| 0 | T | F | median_val = (nums[total_count / 2 - 1] + nums[total_count / 2]) / 2.0; |
| 4883 | 100 | T | F | OUTPUT: |
| 4885 | 50 | T | F | |
| 100 | T | F | |
| 4886 | 50 | T | F | SV* cor(SV* x_sv, SV* y_sv = &PL_sv_undef, const char* method = "pearson") |
| 4893 | 100 | T | F | method); |
| 4894 | 50 | T | F | |
| 4899 | 100 | T | F | AV*restrict x_av = (AV*)SvRV(x_sv); |
| 4902 | 100 | T | F | |
| 4906 | 100 | T | F | SV**restrict fp = av_fetch(x_av, 0, 0); |
| 4908 | 50 | T | F | x_is_matrix = 1; |
| 4910 | 100 | T | F | |
| 100 | T | F | |
| 4911 | 100 | T | F | // --- detect y ---------------------------- |
| 4912 | 50 | T | F | bool has_y = (SvOK(y_sv) && SvROK(y_sv) && |
| 4918 | 50 | T | F | bool y_is_matrix = 0; |
| 4919 | 100 | T | F | if (has_y && ny > 0) { |
| 4920 | 100 | T | F | SV**restrict fp = av_fetch(y_av, 0, 0); |
| 4921 | 100 | T | F | if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV) |
| 4922 | 100 | T | F | y_is_matrix = 1; |
| 4923 | 50 | T | F | } |
| 4936 | 100 | T | F | if (nx < 2) |
| 4947 | 50 | T | F | SV**restrict tv = av_fetch(x_av, i, 0); |
| 50 | T | F | SV**restrict tv = av_fetch(x_av, i, 0); |
| 4948 | 50 | T | F | double val = (tv && SvOK(*tv) && looks_like_number(*tv)) ? SvNV(*tv) : NAN; |
| 4953 | 100 | T | F | } |
| 4955 | 100 | T | F | for (size_t i = 0; i < ny; i++) { |
| 4959 | 100 | T | F | if (!isnan(val)) { |
| 4960 | 100 | T | F | if (isnan(y_first)) y_first = val; |
| 4962 | 100 | T | F | } |
| 4964 | 50 | T | F | |
| 4965 | 100 | T | F | if (x_sd0 || y_sd0) { |
| 4970 | 50 | T | F | |
| 4973 | 50 | T | F | RETVAL = newSVnv(r); |
| 4976 | 100 | T | F | // -- resolve x matrix dimensions |
| 4983 | 100 | T | F | croak("cor: each row of x must be an ARRAY reference"); |
| 4984 | 100 | T | F | |
| 4985 | 100 | T | F | size_t ncols_x = av_len((AV*)SvRV(*xr0)) + 1; |
| 4986 | 100 | T | F | if (ncols_x == 0) croak("cor: x matrix has zero columns"); |
| 4988 | 50 | T | F | size_t nrows = nx; /* observations */ |
| 4992 | 50 | T | F | SV**restrict rv = av_fetch(x_av, i, 0); |
| 5000 | 100 | T | F | SV**restrict rv = av_fetch(y_av, i, 0); |
| 5001 | 100 | T | F | if (!rv || !SvROK(*rv) || SvTYPE(SvRV(*rv)) != SVt_PVAV) |
| 5003 | 100 | T | F | } |
| 5007 | 100 | T | F | double **restrict col_x; |
| 5009 | 100 | T | F | |
| 5015 | 100 | T | F | SV**restrict rv = av_fetch(x_av, i, 0); |
| 5016 | 100 | T | F | AV*restrict row = (AV*)SvRV(*rv); |
| 5019 | 100 | T | F | col_x[j][i] = val; |
| 100 | T | F | col_x[j][i] = val; |
| 5034 | 100 | T | F | double **restrict col_y = NULL; |
| 5037 | 100 | T | F | // 1 = cor(X) — result is symmetric |
| 5039 | 100 | T | F | // cross-correlation: X (nrows × p) vs Y (nrows × q) |
| 100 | T | F | // cross-correlation: X (nrows × p) vs Y (nrows × q) |
| 5042 | 100 | T | F | if (ncols_y == 0) croak("cor: y matrix has zero columns"); |
| 5051 | 50 | T | F | AV*restrict row = (AV*)SvRV(*rv); |
| 5057 | 50 | T | F | else if (val != first) sd0 = 0; |
| 50 | T | F | else if (val != first) sd0 = 0; |
| 5060 | 50 | T | F | if (sd0) { |
| 50 | T | F | if (sd0) { |
| 5063 | 0 | T | F | for (size_t k = 0; k <= j; k++) Safefree(col_y[k]); |
| 5067 | 0 | T | F | } |
| 5071 | 100 | T | F | symmetric = 1; |
| 5075 | 100 | T | F | // -- build cache for symmetric case: compute upper triangle, store results, mirror to lower triangle |
| 5082 | 50 | T | F | rows_out[i] = newAV(); |
| 0 | T | F | rows_out[i] = newAV(); |
| 5102 | 50 | T | F | for (size_t j = 0; j < ncols_x; j++) |
| 5112 | 100 | T | F | } |
| 5113 | 100 | T | F | // push row AVs into result |
| 5114 | 100 | T | F | for (size_t i = 0; i < ncols_x; i++) |
| 5116 | 100 | T | F | Safefree(rows_out); rows_out = NULL; |
| 5121 | 100 | T | F | for (size_t j = 0; j < ncols_y; j++) Safefree(col_y[j]); |
| 5135 | 50 | T | F | size_t data_items = items; |
| 5136 | 0 | T | F | // 1. Parse Options Hash (if it exists as the last argument) |
| 5137 | 0 | T | F | if (items > 0) { |
| 5145 | 100 | T | F | SV*restrict val_sv = *center_sv; |
| 50 | T | F | SV*restrict val_sv = *center_sv; |
| 100 | T | F | SV*restrict val_sv = *center_sv; |
| 50 | T | F | SV*restrict val_sv = *center_sv; |
| 5153 | 50 | T | F | } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) { |
| 5156 | 50 | T | F | do_center_mean = FALSE; center_val = SvNV(val_sv); |
| 5157 | 100 | T | F | } else if (SvTRUE(val_sv)) { |
| 5167 | 100 | T | F | SV*restrict val_sv = *scale_sv; |
| 5174 | 50 | T | F | } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) { |
| 50 | T | F | } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) { |
| 0 | T | F | } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) { |
| 0 | T | F | } else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) { |
| 5180 | 50 | T | F | do_scale_sd = TRUE; |
| 5184 | 100 | T | F | } |
| 5188 | 100 | T | F | // 2. Detect if the input is a Matrix (Array of Arrays) |
| 5189 | 100 | T | F | bool is_matrix = FALSE; |
| 5190 | 50 | T | F | if (data_items == 1) { |
| 5194 | 100 | T | F | if (av_len(av) >= 0) { |
| 5197 | 50 | T | F | is_matrix = TRUE; |
| 5200 | 100 | T | F | } |
| 5207 | 100 | T | F | size_t nrow = av_len(mat_av) + 1, ncol = 0; |
| 50 | T | F | size_t nrow = av_len(mat_av) + 1, ncol = 0; |
| 5212 | 100 | T | F | if (nrow == 0 || ncol == 0) croak("scale requires non-empty matrix"); |
| 5248 | 50 | T | F | croak("scale needs >= 2 rows to calculate standard deviation for a matrix column"); |
| 5253 | 50 | T | F | sum_sq += diff * diff; |
| 5255 | 50 | T | F | col_scale = sqrt(sum_sq / (nrow - 1)); |
| 5257 | 50 | T | F | // Store scaled values back into the new matrix rows |
| 5259 | 50 | T | F | double centered = col_data[r] - col_center; |
| 50 | T | F | double centered = col_data[r] - col_center; |
| 5262 | 50 | T | F | } |
| 5263 | 100 | T | F | Safefree(col_data); |
| 5267 | 0 | T | F | EXTEND(SP, 1); |
| 0 | T | F | EXTEND(SP, 1); |
| 5269 | 0 | T | F | } else { |
| 0 | T | F | } else { |
| 5271 | 0 | T | F | // FLAT LIST MODE: Original functionality |
| 5279 | 0 | T | F | AV*restrict av = (AV*)SvRV(arg); |
| 5282 | 0 | T | F | SV**restrict tv = av_fetch(av, j, 0); |
| 5283 | 0 | T | F | if (tv && SvOK(*tv)) { total_count++; } |
| 5284 | 0 | T | F | } |
| 5286 | 0 | T | F | total_count++; |
| 0 | T | F | total_count++; |
| 0 | T | F | total_count++; |
| 5292 | 0 | T | F | SV*restrict arg = ST(i); |
| 5303 | 100 | T | F | } else if (SvOK(arg)) { |
| 100 | T | F | } else if (SvOK(arg)) { |
| 50 | T | F | } else if (SvOK(arg)) { |
| 5307 | 100 | T | F | } |
| 5308 | 100 | T | F | if (do_center_mean) center_val = sum / total_count; |
| 5309 | 50 | T | F | if (do_scale_sd) { |
| 5317 | 50 | T | F | sum_sq += diff * diff; |
| 5318 | 50 | T | F | } |
| 5319 | 50 | T | F | scale_val = sqrt(sum_sq / (total_count - 1)); |
| 5320 | 50 | T | F | } |
| 0 | T | F | } |
| 5321 | 50 | T | F | EXTEND(SP, total_count); |
| 5322 | 50 | T | F | for (size_t i = 0; i < total_count; i++) { |
| 0 | T | F | for (size_t i = 0; i < total_count; i++) { |
| 5323 | 50 | T | F | double centered = nums[i] - center_val; |
| 0 | T | F | double centered = nums[i] - center_val; |
| 5325 | 50 | T | F | PUSHs(sv_2mortal(newSVnv(final_val))); |
| 5326 | 50 | T | F | } |
| 5328 | 50 | T | F | } |
| 50 | T | F | } |
| 5333 | 100 | T | F | // Basic check: must have an even number of arguments for key => value |
| 5334 | 100 | T | F | if (items % 2 != 0) { |
| 5336 | 100 | T | F | } |
| 5338 | 50 | T | F | size_t nrow = 0, ncol = 0; |
| 50 | T | F | size_t nrow = 0, ncol = 0; |
| 5340 | 100 | T | F | // Parse named arguments |
| 5342 | 50 | T | F | char*restrict key = SvPV_nolen(ST(i)); |
| 5343 | 100 | T | F | SV*restrict val = ST(i + 1); |
| 5353 | 50 | T | F | byrow = SvTRUE(val); |
| 5354 | 100 | T | F | } else { |
| 5370 | 50 | T | F | ncol = 1; |
| 5372 | 50 | T | F | ncol = (data_len + nrow - 1) / nrow; |
| 5374 | 100 | T | F | nrow = (data_len + ncol - 1) / ncol; |
| 5375 | 50 | T | F | } |
| 5380 | 100 | T | F | // Create the matrix (Array of Arrays) |
| 5385 | 50 | T | F | for (r = 0; r < nrow; r++) { |
| 0 | T | F | for (r = 0; r < nrow; r++) { |
| 5386 | 50 | T | F | row_ptrs[r] = newAV(); |
| 0 | T | F | row_ptrs[r] = newAV(); |
| 5394 | 50 | T | F | SV**restrict fetched = av_fetch(data_av, i % data_len, 0); |
| 0 | T | F | SV**restrict fetched = av_fetch(data_av, i % data_len, 0); |
| 5401 | 100 | T | F | c = i / nrow; |
| 5403 | 100 | T | F | av_store(row_ptrs[r], c, val); |
| 5404 | 50 | T | F | } |
| 5406 | 50 | T | F | RETVAL = newRV_noinc((SV*)result_av); |
| 5418 | 100 | T | F | char **restrict terms = NULL, **restrict uniq_terms = NULL, **restrict exp_terms = NULL; |
| 5419 | 100 | T | F | bool *restrict is_dummy = NULL; |
| 5427 | 100 | T | F | HV *restrict data_hoa = NULL; |
| 5437 | 100 | T | F | HE *restrict entry; |
| 5445 | 100 | T | F | else if (strEQ(key, "data")) data_sv = val; |
| 5446 | 100 | T | F | else croak("lm: unknown argument '%s'", key); |
| 5447 | 100 | T | F | } |
| 5450 | 100 | T | F | |
| 50 | T | F | |
| 5454 | 100 | T | F | ref = SvRV(data_sv); |
| 5455 | 100 | T | F | if (SvTYPE(ref) == SVt_PVHV) { |
| 5456 | 50 | T | F | HV *restrict hv = (HV*)ref; |
| 5477 | 100 | T | F | row_hashes[i] = (HV*)SvRV(hv_iterval(hv, entry)); |
| 5481 | 100 | T | F | } |
| 5483 | 50 | T | F | AV *restrict av = (AV*)ref; n = av_len(av) + 1; |
| 5485 | 100 | T | F | Newx(row_hashes, n, HV*); |
| 5486 | 100 | T | F | for (i = 0; i < n; i++) { |
| 5488 | 100 | T | F | if (val && SvROK(*val) && SvTYPE(SvRV(*val)) == SVt_PVHV) { |
| 5489 | 50 | T | F | row_hashes[i] = (HV*)SvRV(*val); |
| 5496 | 50 | T | F | } |
| 5497 | 100 | T | F | } |
| 5498 | 100 | T | F | } else croak("lm: Data must be an Array or Hash reference"); |
| 5499 | 100 | T | F | |
| 5509 | 100 | T | F | for (i = 0; i < n; i++) Safefree(row_names[i]); |
| 5510 | 50 | T | F | Safefree(row_names); if (row_hashes) Safefree(row_hashes); |
| 5527 | 100 | T | F | continue; |
| 5547 | 100 | T | F | continue; |
| 5548 | 50 | T | F | } |
| 5552 | 100 | T | F | } |
| 5554 | 50 | T | F | if (p_idx[0] == '+' && p_idx[1] == '1' && |
| 5557 | 100 | T | F | continue; |
| 5558 | 100 | T | F | } |
| 5560 | 100 | T | F | if (p_idx == rhs) { |
| 5562 | 100 | T | F | if (p_idx[0] == '1' && p_idx[1] == '+') { memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); continue; } |
| 5564 | 50 | T | F | p_idx++; |
| 5565 | 100 | T | F | } |
| 5570 | 50 | T | F | char *restrict p_idx; |
| 5573 | 50 | T | F | if (rhs[0] == '+') memmove(rhs, rhs + 1, strlen(rhs + 1) + 1); |
| 5576 | 100 | T | F | } |
| 5582 | 100 | T | F | while (chunk != NULL) { |
| 5584 | 100 | T | F | AV *cols = get_all_columns(data_hoa, row_hashes, n); |
| 5585 | 100 | T | F | for (size_t c = 0; c <= (size_t)av_len(cols); c++) { |
| 5586 | 100 | T | F | SV **col_sv = av_fetch(cols, c, 0); |
| 5588 | 50 | T | F | const char *col_name = SvPV_nolen(*col_sv); |
| 5594 | 100 | T | F | rhs_len += slen; |
| 5596 | 50 | T | F | } |
| 5598 | 50 | T | F | } |
| 100 | T | F | } |
| 5614 | 100 | T | F | |
| 5615 | 100 | T | F | if (has_intercept) { terms[num_terms++] = savepv("Intercept"); } |
| 5616 | 100 | T | F | |
| 5623 | 100 | T | F | } |
| 5624 | 100 | T | F | char *restrict star = strchr(chunk, '*'); |
| 5627 | 100 | T | F | char *restrict left = chunk; |
| 5631 | 50 | T | F | char *restrict c_r = strchr(right, '^'); |
| 5634 | 100 | T | F | terms[num_terms++] = savepv(right); |
| 5635 | 100 | T | F | size_t inter_len = strlen(left) + strlen(right) + 2; |
| 5639 | 100 | T | F | char *restrict c_chunk = strchr(chunk, '^'); |
| 5644 | 50 | T | F | } |
| 100 | T | F | } |
| 5666 | 100 | T | F | if (is_column_categorical(data_hoa, row_hashes, n, uniq_terms[j])) { |
| 5668 | 50 | T | F | unsigned int num_levels = 0, levels_cap = 8; |
| 50 | T | F | unsigned int num_levels = 0, levels_cap = 8; |
| 5672 | 100 | T | F | if (str_val) { |
| 5674 | 100 | T | F | for (l = 0; l < num_levels; l++) { if (strcmp(levels[l], str_val) == 0) { found = TRUE; break; } } |
| 5676 | 100 | T | F | if (num_levels >= levels_cap) { levels_cap *= 2; Renew(levels, levels_cap, char*); } |
| 5689 | 100 | T | F | Renew(exp_terms, exp_cap, char*); Renew(is_dummy, exp_cap, bool); |
| 5690 | 100 | T | F | Renew(dummy_base, exp_cap, char*); Renew(dummy_level, exp_cap, char*); |
| 5691 | 100 | T | F | } |
| 5693 | 100 | T | F | exp_terms[p_exp] = (char*)safemalloc(t_len); |
| 5700 | 100 | T | F | for (l = 0; l < num_levels; l++) Safefree(levels[l]); |
| 5703 | 50 | T | F | Safefree(levels); |
| 5714 | 100 | T | F | // ======================================================================== |
| 5721 | 100 | T | F | bool row_ok = TRUE; |
| 5722 | 50 | T | F | double *restrict row_x = (double*)safemalloc(p * sizeof(double)); |
| 5725 | 50 | T | F | row_x[j] = 1.0; |
| 50 | T | F | row_x[j] = 1.0; |
| 5727 | 50 | T | F | char *restrict str_val = get_data_string_alloc(data_hoa, row_hashes, i, dummy_base[j]); |
| 5732 | 50 | T | F | } else { |
| 5736 | 100 | T | F | } |
| 5738 | 50 | T | F | |
| 5741 | 50 | T | F | valid_row_names[valid_n] = row_names[i]; |
| 50 | T | F | valid_row_names[valid_n] = row_names[i]; |
| 50 | T | F | valid_row_names[valid_n] = row_names[i]; |
| 50 | T | F | valid_row_names[valid_n] = row_names[i]; |
| 5748 | 50 | T | F | for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms); |
| 50 | T | F | for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms); |
| 5749 | 50 | T | F | for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms); |
| 50 | T | F | for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms); |
| 5750 | 50 | T | F | for (j = 0; j < p_exp; j++) { |
| 50 | T | F | for (j = 0; j < p_exp; j++) { |
| 5751 | 50 | T | F | Safefree(exp_terms[j]); |
| 50 | T | F | Safefree(exp_terms[j]); |
| 5755 | 50 | T | F | Safefree(X); Safefree(Y); Safefree(valid_row_names); |
| 5758 | 50 | T | F | } |
| 5761 | 50 | T | F | // PHASE 5: OLS Math |
| 50 | T | F | // PHASE 5: OLS Math |
| 5764 | 100 | T | F | for (i = 0; i < p; i++) |
| 5765 | 100 | T | F | for (j = 0; j < p; j++) { |
| 5768 | 50 | T | F | XtX[i * p + j] = sum; |
| 50 | T | F | XtX[i * p + j] = sum; |
| 5769 | 50 | T | F | } |
| 50 | T | F | } |
| 5774 | 50 | T | F | XtY[i] = sum; |
| 50 | T | F | XtY[i] = sum; |
| 5779 | 50 | T | F | for (i = 0; i < p; i++) { |
| 5780 | 50 | T | F | if (aliased[i]) { beta[i] = NAN; } |
| 5785 | 50 | T | F | } |
| 5786 | 50 | T | F | } |
| 5787 | 50 | T | F | |
| 50 | T | F | |
| 5788 | 50 | T | F | // ======================================================================== |
| 50 | T | F | // ======================================================================== |
| 5789 | 50 | T | F | // PHASE 6: Metrics & Cleanup |
| 50 | T | F | // PHASE 6: Metrics & Cleanup |
| 5790 | 50 | T | F | // ======================================================================== |
| 50 | T | F | // ======================================================================== |
| 5832 | 50 | T | F | r_squared = 0.0; adj_r_squared = 0.0; |
| 5833 | 100 | T | F | } |
| 5837 | 100 | T | F | av_push(terms_av, newSVpv(exp_terms[j], 0)); |
| 5838 | 100 | T | F | HV *restrict row_hv = newHV(); |
| 5839 | 100 | T | F | if (aliased[j]) { |
| 5840 | 50 | T | F | hv_store(row_hv, "Estimate", 8, newSVpv("NaN", 0), 0); |
| 100 | T | F | hv_store(row_hv, "Estimate", 8, newSVpv("NaN", 0), 0); |
| 5841 | 100 | T | F | hv_store(row_hv, "Std. Error", 10, newSVpv("NaN", 0), 0); |
| 5842 | 100 | T | F | hv_store(row_hv, "t value", 7, newSVpv("NaN", 0), 0); |
| 5843 | 50 | T | F | hv_store(row_hv, "Pr(>|t|)", 8, newSVpv("NaN", 0), 0); |
| 5844 | 0 | T | F | } else { |
| 5845 | 0 | T | F | double se = sqrt(rse_sq * XtX[j * p + j]); |
| 5849 | 100 | T | F | hv_store(row_hv, "Std. Error", 10, newSVnv(se), 0); |
| 50 | T | F | hv_store(row_hv, "Std. Error", 10, newSVnv(se), 0); |
| 5850 | 50 | T | F | hv_store(row_hv, "t value", 7, newSVnv(t_val), 0); |
| 50 | T | F | hv_store(row_hv, "t value", 7, newSVnv(t_val), 0); |
| 5851 | 100 | T | F | hv_store(row_hv, "Pr(>|t|)", 8, newSVnv(p_val), 0); |
| 50 | T | F | hv_store(row_hv, "Pr(>|t|)", 8, newSVnv(p_val), 0); |
| 5852 | 50 | T | F | } |
| 50 | T | F | } |
| 5853 | 100 | T | F | hv_store(summary_hv, exp_terms[j], strlen(exp_terms[j]), newRV_noinc((SV*)row_hv), 0); |
| 50 | T | F | hv_store(summary_hv, exp_terms[j], strlen(exp_terms[j]), newRV_noinc((SV*)row_hv), 0); |
| 5856 | 100 | T | F | hv_store(res_hv, "coefficients", 12, newRV_noinc((SV*)coef_hv), 0); |
| 5857 | 50 | T | F | hv_store(res_hv, "fitted.values", 13, newRV_noinc((SV*)fitted_hv), 0); |
| 5858 | 100 | T | F | hv_store(res_hv, "residuals", 9, newRV_noinc((SV*)resid_hv), 0); |
| 5859 | 50 | T | F | hv_store(res_hv, "df.residual", 11, newSVuv(df_res), 0); |
| 5860 | 50 | T | F | hv_store(res_hv, "rank", 4, newSVuv(final_rank), 0); |
| 5862 | 50 | T | F | hv_store(res_hv, "summary", 7, newRV_noinc((SV*)summary_hv),0); |
| 5866 | 100 | T | F | if (!isnan(f_stat)) { |
| 5867 | 50 | T | F | AV *fstat_av = newAV(); |
| 5868 | 50 | T | F | av_push(fstat_av, newSVnv(f_stat)); |
| 50 | T | F | av_push(fstat_av, newSVnv(f_stat)); |
| 5869 | 100 | T | F | av_push(fstat_av, newSViv(numdf)); |
| 50 | T | F | av_push(fstat_av, newSViv(numdf)); |
| 5870 | 100 | T | F | av_push(fstat_av, newSViv(df_res)); |
| 5871 | 100 | T | F | hv_store(res_hv, "fstatistic", 10, newRV_noinc((SV*)fstat_av), 0); |
| 100 | T | F | hv_store(res_hv, "fstatistic", 10, newRV_noinc((SV*)fstat_av), 0); |
| 5872 | 100 | T | F | hv_store(res_hv, "f.pvalue", 8, newSVnv(f_pvalue), 0); |
| 50 | T | F | hv_store(res_hv, "f.pvalue", 8, newSVnv(f_pvalue), 0); |
| 50 | T | F | hv_store(res_hv, "f.pvalue", 8, newSVnv(f_pvalue), 0); |
| 5873 | 100 | T | F | } |
| 50 | T | F | } |
| 5874 | 100 | T | F | |
| 5876 | 50 | T | F | for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms); |
| 5878 | 50 | T | F | for (j = 0; j < p_exp; j++) { |
| 0 | T | F | for (j = 0; j < p_exp; j++) { |
| 5879 | 100 | T | F | Safefree(exp_terms[j]); |
| 5881 | 100 | T | F | } |
| 5885 | 0 | T | F | if (row_hashes) Safefree(row_hashes); |
| 5887 | 0 | T | F | RETVAL = newRV_noinc((SV*)res_hv); |
| 5889 | 0 | T | F | OUTPUT: |
| 5893 | 0 | T | F | double from |
| 5895 | 0 | T | F | double by |
| 0 | T | F | double by |
| 5896 | 0 | T | F | PPCODE: |
| 5898 | 0 | T | F | //Handle the zero 'by' case |
| 5902 | 0 | T | F | mPUSHn(from); |
| 5904 | 0 | T | F | } else { |
| 5906 | 0 | T | F | } |
| 5918 | 100 | T | F | size_t n_elements = (n_elements_d + 1e-10) + 1; |
| 100 | T | F | size_t n_elements = (n_elements_d + 1e-10) + 1; |
| 5920 | 100 | T | F | EXTEND(SP, n_elements); |
| 100 | T | F | EXTEND(SP, n_elements); |
| 5921 | 100 | T | F | for (size_t i = 0; i < n_elements; i++) { |
| 5935 | 50 | T | F | int arg_start = 0; |
| 100 | T | F | int arg_start = 0; |
| 5937 | 100 | T | F | // Check if the first argument is a simple integer (rnorm(33)) |
| 5939 | 50 | T | F | n = (unsigned int)SvUV(ST(0)); |
| 5943 | 100 | T | F | // --- Parse remaining named arguments from the flat stack --- |
| 50 | T | F | // --- Parse remaining named arguments from the flat stack --- |
| 5944 | 100 | T | F | if ((items - arg_start) % 2 != 0) { |
| 5945 | 50 | T | F | croak("Usage: rnorm(n), rnorm(n => 10, mean => 0, sd => 1), or rnorm(33, mean => 0)"); |
| 5949 | 100 | T | F | const char* restrict key = SvPV_nolen(ST(i)); |
| 5952 | 100 | T | F | if (strEQ(key, "n")) n = (unsigned int)SvUV(val); |
| 5953 | 50 | T | F | else if (strEQ(key, "mean")) mean = SvNV(val); |
| 5954 | 0 | T | F | else if (strEQ(key, "sd")) sd = SvNV(val); |
| 5958 | 100 | T | F | if (sd < 0.0) croak("rnorm: standard deviation must be non-negative"); |
| 50 | T | F | if (sd < 0.0) croak("rnorm: standard deviation must be non-negative"); |
| 50 | T | F | if (sd < 0.0) croak("rnorm: standard deviation must be non-negative"); |
| 5973 | 100 | T | F | double mul = sqrt(-2.0 * log(s) / s); |
| 5974 | 50 | T | F | // Box-Muller generates two independent values per iteration |
| 50 | T | F | // Box-Muller generates two independent values per iteration |
| 5982 | 100 | T | F | } |
| 5984 | 50 | T | F | RETVAL |
| 50 | T | F | RETVAL |
| 5988 | 50 | T | F | SV* formula_sv |
| 5992 | 50 | T | F | SV *restrict orig_data_sv = data_sv; |
| 5998 | 100 | T | F | if (!formula_sv || !SvOK(formula_sv) || SvCUR(formula_sv) == 0) { |
| 6000 | 50 | T | F | croak("aov: Without a formula, data must be a HashRef of ArrayRefs (mimicking R's named list)"); |
| 6005 | 100 | T | F | HV *restrict stacked_hv = newHV(); |
| 6007 | 50 | T | F | AV *restrict grp_av = newAV(); |
| 50 | T | F | AV *restrict grp_av = newAV(); |
| 50 | T | F | AV *restrict grp_av = newAV(); |
| 6021 | 50 | T | F | av_push(val_av, newSVsv(*v)); |
| 50 | T | F | av_push(val_av, newSVsv(*v)); |
| 50 | T | F | av_push(val_av, newSVsv(*v)); |
| 6023 | 50 | T | F | } |
| 50 | T | F | } |
| 50 | T | F | } |
| 6030 | 50 | T | F | |
| 6031 | 50 | T | F | hv_stores(stacked_hv, "Value", newRV_noinc((SV*)val_av)); |
| 6040 | 100 | T | F | |
| 6043 | 50 | T | F | |
| 50 | T | F | |
| 50 | T | F | |
| 6044 | 50 | T | F | char **restrict terms = NULL, **restrict uniq_terms = NULL, **restrict exp_terms = NULL, **restrict parent_term = NULL; |
| 50 | T | F | char **restrict terms = NULL, **restrict uniq_terms = NULL, **restrict exp_terms = NULL, **restrict parent_term = NULL; |
| 6049 | 100 | T | F | size_t n = 0, valid_n = 0, i, j; |
| 6068 | 50 | T | F | HV*restrict hv = (HV*)ref; |
| 50 | T | F | HV*restrict hv = (HV*)ref; |
| 6070 | 0 | T | F | entry = hv_iternext(hv); |
| 6071 | 0 | T | F | if (entry) { |
| 0 | T | F | if (entry) { |
| 6074 | 0 | T | F | data_hoa = hv; |
| 6087 | 100 | T | F | row_names[i] = savepv(hv_iterkey(entry, &len)); |
| 6096 | 100 | T | F | Newx(row_names, n, char*); |
| 6097 | 50 | T | F | Newx(row_hashes, n, HV*); |
| 6104 | 50 | T | F | row_names[i] = savepv(buf); |
| 6124 | 100 | T | F | croak("aov: invalid formula, missing '~'"); |
| 6125 | 50 | T | F | } |
| 50 | T | F | } |
| 6132 | 50 | T | F | while ((p_idx = strstr(rhs, "+0")) != NULL) { has_intercept = FALSE; memmove(p_idx, p_idx + 2, strlen(p_idx + 2) + 1); } |
| 6162 | 50 | T | F | } |
| 50 | T | F | } |
| 50 | T | F | } |
| 6168 | 50 | T | F | if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; } |
| 50 | T | F | if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; } |
| 50 | T | F | if (rhs_len > 0) { strcat(rhs_expanded, "+"); rhs_len++; } |
| 6174 | 50 | T | F | } |
| 6179 | 100 | T | F | Newx(is_dummy, exp_cap, bool); Newx(is_interact, exp_cap, bool); |
| 6183 | 50 | T | F | if (strlen(rhs_expanded) > 0) { |
| 6184 | 50 | T | F | chunk = strtok(rhs_expanded, "+"); |
| 6185 | 100 | T | F | while (chunk != NULL) { |
| 6186 | 50 | T | F | if (num_terms >= term_cap - 3) { |
| 0 | T | F | if (num_terms >= term_cap - 3) { |
| 6187 | 0 | T | F | term_cap *= 2; |
| 6192 | 50 | T | F | *star = '\0'; |
| 50 | T | F | *star = '\0'; |
| 50 | T | F | *star = '\0'; |
| 6194 | 50 | T | F | char *right = star + 1; |
| 50 | T | F | char *right = star + 1; |
| 50 | T | F | char *right = star + 1; |
| 6197 | 50 | T | F | char *restrict c_r = strchr(right, '^'); if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0'; |
| 50 | T | F | char *restrict c_r = strchr(right, '^'); if (c_r && strncmp(right, "I(", 2) != 0) *c_r = '\0'; |
| 6199 | 50 | T | F | terms[num_terms++] = savepv(right); |
| 50 | T | F | terms[num_terms++] = savepv(right); |
| 50 | T | F | terms[num_terms++] = savepv(right); |
| 6210 | 100 | T | F | } |
| 6212 | 50 | T | F | for (i = 0; i < num_terms; i++) { |
| 50 | T | F | for (i = 0; i < num_terms; i++) { |
| 50 | T | F | for (i = 0; i < num_terms; i++) { |
| 6214 | 50 | T | F | for (size_t k = 0; k < num_uniq; k++) { |
| 50 | T | F | for (size_t k = 0; k < num_uniq; k++) { |
| 6225 | 100 | T | F | if (p_exp + 64 >= exp_cap) { |
| 6227 | 50 | T | F | Renew(exp_terms, exp_cap, char*); Renew(parent_term, exp_cap, char*); |
| 50 | T | F | Renew(exp_terms, exp_cap, char*); Renew(parent_term, exp_cap, char*); |
| 50 | T | F | Renew(exp_terms, exp_cap, char*); Renew(parent_term, exp_cap, char*); |
| 6229 | 50 | T | F | Renew(dummy_base, exp_cap, char*); Renew(dummy_level, exp_cap, char*); |
| 50 | T | F | Renew(dummy_base, exp_cap, char*); Renew(dummy_level, exp_cap, char*); |
| 6238 | 100 | T | F | p_exp++; |
| 6239 | 100 | T | F | continue; |
| 6246 | 100 | T | F | left[colon - uniq_terms[j]] = '\0'; |
| 6255 | 50 | T | F | |
| 6257 | 50 | T | F | Safefree(l_indices); Safefree(r_indices); |
| 6264 | 50 | T | F | Renew(exp_terms, exp_cap, char*); Renew(parent_term, exp_cap, char*); |