Stats-LikeR
view release on metacpan or search on metacpan
SV *restrict tv = src[j];
if (tv && SvOK(tv)) moment_push(acc, SvNV(tv));
else croak("%s: undefined value at array ref index %" UVuf
" (argument %" UVuf ")", fname, (UV)j, (UV)argi);
}
}
/*Walk an argument list of numbers, array refs of numbers and 'type'/'x'
named pairs. Shared so that skew() and kurtosis() cannot drift apart on
what they accept. A named key is recognised only when the SV is a string
that is not a number, so it can never swallow a data value -- and anything
else that looks like a bareword is a typo worth reporting rather than
silently averaging in as zero.*/
static void moment_args(pTHX_ SV **restrict args, size_t items,
const char *restrict fname,
moment_acc *restrict acc, IV *restrict type) {
for (size_t i = 0; i < items; i++) {
SV *restrict arg = args[i];
if (arg && SvPOK(arg) && !SvROK(arg) && !looks_like_number(arg)) {
const char *restrict key = SvPV_nolen(arg);
const bool is_type = strEQ(key, "type");
if (!is_type && !strEQ(key, "x"))
croak("%s: unknown argument '%s' (expected numbers, array "
"refs, x => \\@data or type => 1|2|3)", fname, key);
if (i + 1 >= items)
croak("%s: '%s' needs a value", fname, key);
SV *restrict val = args[++i];
if (is_type) {
if (!SvOK(val) || !looks_like_number(val))
croak("%s: type must be 1, 2 or 3", fname);
*type = SvIV(val);
if (*type < 1 || *type > 3)
croak("%s: type must be 1, 2 or 3, not %" IVdf, fname, *type);
} else {
if (!SvROK(val) || SvTYPE(SvRV(val)) != SVt_PVAV)
croak("%s: 'x' must be an array reference", fname);
moment_av(aTHX_ (AV*)SvRV(val), i, fname, acc);
}
} else if (arg && SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
moment_av(aTHX_ (AV*)SvRV(arg), i, fname, acc);
} else if (arg && SvOK(arg)) {
moment_push(acc, SvNV(arg));
} else {
croak("%s: undefined value at argument index %" UVuf, fname, (UV)i);
}
}
}
// --- XS SECTION ---
MODULE = Stats::LikeR PACKAGE = Stats::LikeR
void
_interp_column_xs(vals_ref, x_ref, method, order_sv, dir, limit_sv, area_sv)
SV *vals_ref
SV *x_ref
const char *method
SV *order_sv
const char *dir
SV *limit_sv
SV *area_sv
PPCODE:
if (!(SvROK(vals_ref) && SvTYPE(SvRV(vals_ref)) == SVt_PVAV))
croak("_interp_column_xs: values must be an array reference");
if (!(SvROK(x_ref) && SvTYPE(SvRV(x_ref)) == SVt_PVAV))
croak("_interp_column_xs: x must be an array reference");
ENTER; SAVETMPS;
ip_fill_column(aTHX_ (AV *)SvRV(vals_ref), (AV *)SvRV(x_ref),
method, order_sv, dir, limit_sv, area_sv);
FREETMPS; LEAVE;
XSRETURN_EMPTY;
SV *_cols_select(df, shape, spec)
SV *df
IV shape
SV *spec
PREINIT:
SV *restrict retval; AV *restrict spec_av; SSize_t n, i;
CODE:
{
spec_av = (AV *)SvRV(spec);
n = av_len(spec_av) + 1;
if (shape == 3) { // ---- AoA ----
IV *restrict idx; Newx(idx, n > 0 ? n : 1, IV);
for (i = 0; i < n; i++) { SV **e = av_fetch(spec_av, i, 0); idx[i] = SvIV(*e); }
AV *restrict src = (AV *)SvRV(df); SSize_t R = av_len(src) + 1;
AV *restrict out = newAV(); if (R > 0) av_extend(out, R - 1);
for (i = 0; i < R; i++) {
SV **restrict rp = av_fetch(src, i, 0); AV *inner;
if (rp && *rp && SvROK(*rp) && SvTYPE(SvRV(*rp)) == SVt_PVAV)
inner = rowA_select(aTHX_ (AV *)SvRV(*rp), idx, n);
else
inner = rowA_select(aTHX_ NULL, idx, n);
av_store(out, i, newRV_noinc((SV *)inner));
}
Safefree(idx);
retval = sv_2mortal(newRV_noinc((SV *)out));
} else {
SV **restrict keys; Newx(keys, n > 0 ? n : 1, SV *);
for (i = 0; i < n; i++) { SV **e = av_fetch(spec_av, i, 0); keys[i] = *e; }
if (shape == 1) { // ---- AoH ----
AV *restrict src = (AV *)SvRV(df); SSize_t R = av_len(src) + 1;
AV *restrict out = newAV(); if (R > 0) av_extend(out, R - 1);
for (i = 0; i < R; i++) {
SV **restrict rp = av_fetch(src, i, 0); HV *inner;
if (rp && *rp && SvROK(*rp) && SvTYPE(SvRV(*rp)) == SVt_PVHV)
inner = row_select(aTHX_ (HV *)SvRV(*rp), keys, n);
else
inner = row_select(aTHX_ NULL, keys, n);
av_store(out, i, newRV_noinc((SV *)inner));
}
retval = sv_2mortal(newRV_noinc((SV *)out));
} else { // ---- HoH ----
HV *restrict src = (HV *)SvRV(df); HV *out = newHV();
hv_iterinit(src); HE *restrict he;
while ((he = hv_iternext(src))) {
STRLEN kl; char *restrict kp = HePV(he, kl); I32 sk = HeUTF8(he) ? -(I32)kl : (I32)kl;
SV *restrict rv = HeVAL(he); HV *inner;
if (rv && SvROK(rv) && SvTYPE(SvRV(rv)) == SVt_PVHV)
inner = row_select(aTHX_ (HV *)SvRV(rv), keys, n);
else
inner = row_select(aTHX_ NULL, keys, n);
if (shape == 4) { // ---- HoA ----
HV *restrict src = (HV *)SvRV(df);
HV *restrict out = newHV();
HE *restrict he; hv_iterinit(src);
while ((he = hv_iternext(src))) {
SV *restrict v = HeVAL(he);
AV *restrict col = (SvROK(v) && SvTYPE(SvRV(v)) == SVt_PVAV) ? (AV *)SvRV(v) : NULL;
const SSize_t cfill = col ? AvFILLp(col) : -1;
SV **restrict ca = col ? AvARRAY(col) : NULL;
AV *restrict nc = newAV();
if (nsurv > 0) {
av_extend(nc, nsurv - 1);
SV **restrict na = AvARRAY(nc);
for (SSize_t k = 0; k < nsurv; k++) {
SV *restrict c = surv[k] <= cfill ? ca[surv[k]] : NULL;
na[k] = c ? newSVsv(c) : newSV(0);
}
AvFILLp(nc) = nsurv - 1;
}
STRLEN kl; char *kp = HePV(he, kl); I32 sk = HeUTF8(he) ? -(I32)kl : (I32)kl;
(void)hv_store(out, kp, sk, newRV_noinc((SV *)nc), HeHASH(he));
}
retval = sv_2mortal(newRV_noinc((SV *)out));
} else { // AoA / AoH
AV *restrict src = (AV *)SvRV(df);
const SSize_t sfill = AvFILLp(src);
SV **restrict sa = AvARRAY(src);
AV *restrict out = newAV();
if (nsurv > 0) {
av_extend(out, nsurv - 1);
SV **restrict oa = AvARRAY(out);
for (SSize_t k = 0; k < nsurv; k++) {
SV *restrict rv = surv[k] <= sfill ? sa[surv[k]] : NULL;
oa[k] = newSVsv(rv ? rv : &PL_sv_undef);
}
AvFILLp(out) = nsurv - 1;
}
retval = sv_2mortal(newRV_noinc((SV *)out));
}
RETVAL = SvREFCNT_inc(retval);
LEAVE; //dd_ctx_free releases the rest
}
OUTPUT:
RETVAL
void anova(...)
PROTOTYPE: $@
PREINIT:
SV *restrict data;
char *lhs = NULL, *rhs = NULL;
HV *restrict hoa = NULL, *restrict result = NULL;
HV **restrict rows = NULL;
AnTerm *terms = NULL;
AnFac *facs = NULL;
size_t nterms = 0, tcap = 0, nfac = 0, fcap = 0;
bool *restrict complete = NULL, *restrict aliased = NULL;
size_t *restrict ridx = NULL, *rank_map = NULL;
size_t n = 0, n_used = 0, p, rank;
NV **restrict X = NULL, *restrict y = NULL, rss, msres;
IV dfres;
PPCODE:
{
if (items < 2)
croak("anova: usage anova(\\%%data, 'response ~ terms' [, 'model2', ...])");
data = ST(0);
if (items > 2) {
/*nested model comparison *
anova(\%data, 'y ~ a', 'y ~ a + b', ...) -> ArrayRef table.*/
size_t nform = (size_t)items - 1;
char **restrict lhss = NULL, **rhss = NULL;
Newxz(lhss, nform, char*);
Newxz(rhss, nform, char*);
//---- parse every formula
for (size_t fi = 0; fi < nform; fi++) {
SV *restrict fsv = ST(1 + fi);
if (!(SvPOK(fsv) || SvOK(fsv))) {
anova_free_formulas(aTHX_ lhss, rhss, nform);
croak("anova: model argument %" UVuf " must be a formula string", (UV)(fi + 1));
}
if (!parse_formula(SvPV_nolen(fsv), &lhss[fi], &rhss[fi])) {
anova_free_formulas(aTHX_ lhss, rhss, nform);
croak("anova: could not parse formula %" UVuf " (need 'response ~ terms')", (UV)(fi + 1));
}
}
// ---- resolve data form + row count (response 1 length)
if (!SvROK(data)) {
anova_free_formulas(aTHX_ lhss, rhss, nform);
croak("anova: first argument must be a hash or array reference");
}
{
SV *rv = SvRV(data);
if (SvTYPE(rv) == SVt_PVHV) {
hoa = (HV*)rv;
SV **restrict col = hv_fetch(hoa, lhss[0], (I32)strlen(lhss[0]), 0);
if (col && SvROK(*col) && SvTYPE(SvRV(*col)) == SVt_PVAV)
n = (size_t)(av_len((AV*)SvRV(*col)) + 1);
else {
hv_iterinit(hoa);
HE *restrict e;
while ((e = hv_iternext(hoa))) {
SV *v = hv_iterval(hoa, e);
if (SvROK(v) && SvTYPE(SvRV(v)) == SVt_PVAV) {
size_t l = (size_t)(av_len((AV*)SvRV(v)) + 1);
if (l > n) n = l;
}
}
}
} else if (SvTYPE(rv) == SVt_PVAV) {
AV *top = (AV*)rv;
n = (size_t)(av_len(top) + 1);
Newx(rows, n ? n : 1, HV*);
for (size_t i = 0; i < n; i++) {
SV **ep = av_fetch(top, i, 0);
if (!(ep && SvROK(*ep) && SvTYPE(SvRV(*ep)) == SVt_PVHV)) {
Safefree(rows);
anova_free_formulas(aTHX_ lhss, rhss, nform);
croak("anova: element %" UVuf " is not a hash reference", (UV)i);
}
rows[i] = (HV*)SvRV(*ep);
X[r][terms[t].start + c] = v;
}
}
}
// sequential QR (X, y overwritten in place)
Newx(aliased, p, bool);
Newx(rank_map, p, size_t);
for (size_t k = 0; k < p; k++) rank_map[k] = 0;
apply_householder_aov(X, y, n_used, p, aliased, rank_map);
rank = 0;
for (size_t k = 0; k < p; k++) if (!aliased[k]) rank++;
rss = 0.0;
for (size_t r = rank; r < n_used; r++) rss += y[r] * y[r];
dfres = (IV)n_used - (IV)rank;
msres = dfres > 0 ? rss / (NV)dfres : NAN;
// assemble term-keyed table
result = newHV();
for (size_t t = 0; t < nterms; t++) {
NV ss = 0.0; IV df = 0;
for (size_t k = terms[t].start; k < terms[t].start + terms[t].width; k++)
if (!aliased[k]) { ss += y[rank_map[k]] * y[rank_map[k]]; df++; }
HV *restrict in = newHV();
(void)hv_store(in, "Df", 2, newSViv(df), 0);
(void)hv_store(in, "Sum Sq", 6, newSVnv(ss), 0);
if (df > 0) {
(void)hv_store(in, "Mean Sq", 7, newSVnv(ss / (NV)df), 0);
if (dfres > 0 && rss > 0.0) {
NV F = (ss / (NV)df) / msres;
(void)hv_store(in, "F value", 7, newSVnv(F), 0);
(void)hv_store(in, "Pr(>F)", 6, newSVnv(pf_upper(F, (NV)df, (NV)dfres)), 0);
}
}
(void)hv_store(result, terms[t].name, (I32)strlen(terms[t].name),
newRV_noinc((SV*)in), 0);
}
{
HV *restrict in = newHV();
(void)hv_store(in, "Df", 2, newSViv(dfres), 0);
(void)hv_store(in, "Sum Sq", 6, newSVnv(rss), 0);
if (dfres > 0) (void)hv_store(in, "Mean Sq", 7, newSVnv(msres), 0);
(void)hv_store(result, "Residuals", 9, newRV_noinc((SV*)in), 0);
}
// teardown
for (size_t r = 0; r < n_used; r++) Safefree(X[r]);
Safefree(X); Safefree(y);
Safefree(aliased); Safefree(rank_map); Safefree(ridx); Safefree(complete);
anova_free_terms(aTHX_ terms, nterms);
anova_free_facs(aTHX_ facs, nfac);
Safefree(rows);
safefree(lhs); safefree(rhs);
XPUSHs(sv_2mortal(newRV_noinc((SV*)result)));
}
}
void rank(...)
PROTOTYPE: @
PPCODE:
int ties = RANK_AVERAGE;
int nalast = NALAST_TRUE;
/* ---- locate trailing "key => value" options -------------
Options begin at the first plain-string arg equal to a
known option name; everything before it is data.*/
int opt_start = items;
for (int i = 0; i < items; i++) {
SV *a = ST(i);
if (SvOK(a) && !SvROK(a) && SvPOK(a)) {
STRLEN klen;
const char *k = SvPV_const(a, klen);
if ((klen == 11 && strEQ(k, "ties.method")) ||
(klen == 7 && strEQ(k, "na.last"))) {
opt_start = i;
break;
}
}
}
if (((items - opt_start) & 1) != 0)
croak("rank: named options must be key => value pairs");
for (int i = opt_start; i < items; i += 2) {
STRLEN klen, vlen;
const char *k = SvPV_const(ST(i), klen);
SV *vsv = ST(i + 1);
if (strEQ(k, "ties.method")) {
if (!SvOK(vsv))
croak("rank: ties.method cannot be undef");
const char *v = SvPV_const(vsv, vlen);
if (strEQ(v, "average")) ties = RANK_AVERAGE;
else if (strEQ(v, "first")) ties = RANK_FIRST;
else if (strEQ(v, "last")) ties = RANK_LAST;
else if (strEQ(v, "random")) ties = RANK_RANDOM;
else if (strEQ(v, "max")) ties = RANK_MAX;
else if (strEQ(v, "min")) ties = RANK_MIN;
else croak("rank: unknown ties.method '%s' "
"(average, first, last, random, max, min)", v);
} else if (strEQ(k, "na.last")) {
if (!SvOK(vsv)) {
nalast = NALAST_DROP; // undef => R's NA
} else {
const char *v = SvPV_const(vsv, vlen);
if (strEQ(v, "keep")) nalast = NALAST_KEEP;
else if (strEQ(v, "na") || strEQ(v, "NA")) nalast = NALAST_DROP;
else if (strEQ(v, "false") || strEQ(v, "FALSE")
|| strEQ(v, "F") || strEQ(v, "0")) nalast = NALAST_FALSE;
else if (strEQ(v, "true") || strEQ(v, "TRUE")
|| strEQ(v, "T") || strEQ(v, "1")) nalast = NALAST_TRUE;
else croak("rank: unknown na.last '%s' "
"(true, false, keep, na)", v);
}
} else {
croak("rank: unknown option '%s' (ties.method, na.last)", k);
}
}
// ---- count total data elements --------------------------
size_t N = 0;
PVAL = bt_pbinom_lower(y - 1, n, p) + bt_pbinom_upper(x - 1, n, p);
}
}
}
if (PVAL > 1.0) PVAL = 1.0;
// confidence interval (Clopper-Pearson)
NV ci_lo, ci_hi;
if (strEQ(alternative, "less")) {
ci_lo = 0.0;
ci_hi = bt_pU(1.0 - conf_level, x, n);
} else if (strEQ(alternative, "greater")) {
ci_lo = bt_pL(1.0 - conf_level, x, n);
ci_hi = 1.0;
} else {
NV a = (1.0 - conf_level) / 2.0;
ci_lo = bt_pL(a, x, n);
ci_hi = bt_pU(a, x, n);
}
// ---- htest-style result ----
HV *restrict ret = newHV();
hv_stores(ret, "method", newSVpv("Exact binomial test", 0));
hv_stores(ret, "alternative", newSVpv(alternative, 0));
hv_stores(ret, "statistic", newSViv(x)); //number of successes
hv_stores(ret, "parameter", newSViv(n)); //number of trials
hv_stores(ret, "estimate", newSVnv((NV)x / (NV)n)); //probability of success
hv_stores(ret, "null_value", newSVnv(p));
hv_stores(ret, "p_value", newSVnv(PVAL));
hv_stores(ret, "conf_level", newSVnv(conf_level));
AV *restrict ci = newAV();
av_push(ci, newSVnv(ci_lo));
av_push(ci, newSVnv(ci_hi));
hv_stores(ret, "conf_int", newRV_noinc((SV *)ci));
RETVAL = newRV_noinc((SV *)ret);
}
OUTPUT:
RETVAL
BOOT:
newXS("Stats::LikeR::__cs_uninit_catcher", cs_uninit_catcher, __FILE__);
void csort(...)
PREINIT:
SV *restrict data = NULL, *restrict by = NULL, *restrict output = NULL;
cs_shape in_shape = CS_AOH, out_shape = CS_AOH;
bool is_hoh = 0, is_code = 0;
const char *restrict colname = NULL;
STRLEN collen = 0;
IV aoa_col = 0; // AoA: parsed non-negative column index
const char *restrict rowname_col = NULL; // HoH: row-name column name
STRLEN rowname_len = 0;
CV *restrict cmp_cv = NULL;
AV *restrict src_av = NULL; // AoH / AoA input
HV *restrict src_hv = NULL; // HoA / HoH input
SSize_t n = 0;
size_t *restrict idx = NULL, *tmp = NULL;
SV **restrict rowrefs = NULL; // coderef mode: row ref per index
SV **restrict colkeys = NULL; // HoA: column key SVs
AV **restrict colavs = NULL; // HoA: column AVs
size_t ncols = 0;
SV *restrict result = NULL;
PPCODE:
{
// ---- own the usage message (variadic: xsubpp won't invent one)
if (items < 2 || items > 4)
croak("Usage: csort($df, 'column.name', 'HoA')\n"
" or csort($df, sub { $b->{'No.'} <=> $a->{'No.'} }, 'hoa')\n"
" or csort($aoa, 0, 'aoa') # array-of-arrays, integer column\n"
" (optional 4th arg names the row-name column when sorting a "
"HoH; default 'row.name')");
data = ST(0);
by = ST(1);
output = (items >= 3) ? ST(2) : &PL_sv_undef;
if (items >= 4 && SvOK(ST(3)))
rowname_col = SvPV(ST(3), rowname_len);
else {
rowname_col = "row.name";
rowname_len = 8;
}
ENTER; // scope for SAVEFREEPV / SAVESPTR cleanups
SAVETMPS; // reap transient synthesized rows and mortals here
// classify $by: coderef comparator vs column name/index
if (SvROK(by) && SvTYPE(SvRV(by)) == SVt_PVCV) {
is_code = 1;
cmp_cv = (CV *)SvRV(by);
} else if (SvOK(by) && !SvROK(by)) {
is_code = 0;
colname = SvPV(by, collen);
} else {
croak("csort: second argument must be a column name (e.g. 'No.'), an "
"integer column index for an AoA, or a comparator code-ref "
"using $a and $b, e.g. sub { $b->{'No.'} <=> $a->{'No.'} }");
}
//---- classify $data: AoH/AoA (arrayref) vs HoA/HoH (hashref) ------
if (!SvROK(data))
croak("csort: first argument must be an array-ref (AoH or AoA) or "
"hash-ref (HoA or HoH); Usage: csort($df, 'column.name', 'HoA')");
if (SvTYPE(SvRV(data)) == SVt_PVAV) {
src_av = (AV *)SvRV(data);
n = av_len(src_av) + 1;
in_shape = CS_AOH; //default; refine by peeking at row 0
if (n > 0) {
SV **restrict rp = av_fetch(src_av, 0, 0);
if (rp && *rp && SvROK(*rp)
&& SvTYPE(SvRV(*rp)) == SVt_PVAV)
in_shape = CS_AOA; //first row is an arrayref => AoA
}
} else if (SvTYPE(SvRV(data)) == SVt_PVHV) {
src_hv = (HV *)SvRV(data);
hv_iterinit(src_hv);
HE *restrict he = hv_iternext(src_hv);
if (!he) {
in_shape = CS_HOA; //empty hash defaults to HoA path
} else {
SV *restrict val = HeVAL(he);
if (SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV)
is_hoh = 1;
else
in_shape = CS_HOA;
}
} else {
SV *restrict rv = hv_iterval(in_hv, e);
if (!SvROK(rv) || SvTYPE(SvRV(rv)) != SVt_PVHV) croak("hoh2hoa: every value must be a hash ref (hash of hashes)");
av_push(rows_av, newSVsv(hv_iterkeysv(e)));
}
}
SSize_t nrows = av_len(rows_av) + 1;
if (nrows > 1) qsort(AvARRAY(rows_av), (size_t)nrows, sizeof(SV*), h2h_keycmp);
/* 4. discover the union of inner keys. Each new column gets an empty array
in the result straight away so step 5 can just push into it.*/
{
HE *restrict e;
hv_iterinit(in_hv);
while ((e = hv_iternext(in_hv))) {
HV *restrict row = (HV*)SvRV(hv_iterval(in_hv, e));
HE *restrict ie;
hv_iterinit(row);
while ((ie = hv_iternext(row))) {
SV *restrict ck = hv_iterkeysv(ie);
if (!hv_exists_ent(seen, ck, 0)) {
(void)hv_store_ent(seen, ck, &PL_sv_yes, 0);
av_push(cols_av, newSVsv(ck));
(void)hv_store_ent(out_hv, ck, newRV_noinc((SV*)newAV()), 0);
}
}
}
}
SSize_t ncols = av_len(cols_av) + 1;
/* 5. walk the rows in sorted order; for every column push the cell (a copy)
or the fill value, so each column ends up exactly nrows long.*/
for (SSize_t r = 0; r < nrows; r++) {
SV *restrict rk = *av_fetch(rows_av, r, 0);
HE *restrict rhe = hv_fetch_ent(in_hv, rk, 0, 0);
HV *restrict row = (HV*)SvRV(HeVAL(rhe));
for (SSize_t c = 0; c < ncols; c++) {
SV *restrict ck = *av_fetch(cols_av, c, 0);
HE *restrict che = hv_fetch_ent(row, ck, 0, 0);
SV *restrict src = che ? HeVAL(che) : NULL;
SV *restrict cell = (src && SvOK(src)) ? newSVsv(src) : (fill ? newSVsv(fill) : newSV(0));
HE *restrict colhe = hv_fetch_ent(out_hv, ck, 0, 0);
av_push((AV*)SvRV(HeVAL(colhe)), cell);
}
}
// 6. optional row-names column: the sorted labels under the requested name.
if (rn_sv) {
if (hv_exists_ent(out_hv, rn_sv, 0)) croak("hoh2hoa: row.names column '%s' collides with an existing column", SvPV_nolen(rn_sv));
AV *restrict rn_av = newAV();
for (SSize_t r = 0; r < nrows; r++) av_push(rn_av, newSVsv(*av_fetch(rows_av, r, 0)));
(void)hv_store_ent(out_hv, rn_sv, newRV_noinc((SV*)rn_av), 0);
}
// 7. tidy up the scratch structures (the result keeps its own copies).
SvREFCNT_dec((SV*)rows_av);
SvREFCNT_dec((SV*)cols_av);
SvREFCNT_dec((SV*)seen);
RETVAL = newRV_noinc((SV*)out_hv);
}
OUTPUT:
RETVAL
void filter(...)
PPCODE:
{
if (items < 2)
croak("Usage: filter($df, $code [, 'output.type' => 'aoh'|'hoa'])");
SV *restrict df = ST(0);
SV *restrict predarg = ST(1);
const char *restrict otype = NULL;
if (items == 3) {
otype = SvPV_nolen(ST(2));
} else if (items == 4) {
const char *restrict key = SvPV_nolen(ST(2));
if (strNE(key, "output.type") && strNE(key, "out") && strNE(key, "output_type"))
croak("filter: unknown option '%s' (expected 'output.type')", key);
otype = SvPV_nolen(ST(3));
} else if (items > 4) {
croak("Usage: filter($df, $code [, 'output.type' => 'aoh'|'hoa'])");
}
int want = 0; // 0 = preserve input shape
if (otype) {
if (strEQ(otype, "aoh")) want = FLT_AOH;
else if (strEQ(otype, "hoa")) want = FLT_HOA;
else croak("filter: output.type must be 'aoh' or 'hoa' (got '%s')", otype);
}
if (!df || !SvROK(df))
croak("filter: first argument must be a data frame (AoH, HoA, or HoH reference)");
/*The predicate is a CODE ref, or a col() object carrying a CODE ref in
its {code} field; either way we end up calling a single CV per row --
unless the col() object also carries a {plan}, which is the same test as
data and runs in C without touching perl at all.*/
SV *restrict code = NULL;
SV *restrict plan = NULL;
if (predarg && SvROK(predarg) && SvTYPE(SvRV(predarg)) == SVt_PVCV) {
code = predarg;
} else if (predarg && sv_isobject(predarg)
&& sv_derived_from(predarg, "Stats::LikeR::col")) {
SV **restrict cp = hv_fetchs((HV*)SvRV(predarg), "code", 0);
if (!cp || !*cp || !SvROK(*cp) || SvTYPE(SvRV(*cp)) != SVt_PVCV)
croak("filter: incomplete col() predicate -- a bare column needs a comparison, e.g. col('x') > 0");
code = *cp;
SV **restrict pp = hv_fetchs((HV*)SvRV(predarg), "plan", 0);
if (pp && *pp && SvROK(*pp)) plan = *pp;
} else {
croak("filter: predicate must be a CODE reference or a col() expression");
}
SV *restrict ref = SvRV(df);
int in_shape;
HV *restrict inhv = NULL;
AV *restrict inav = NULL;
if (SvTYPE(ref) == SVt_PVAV) {
in_shape = FLT_AOH; inav = (AV*)ref;
} else if (SvTYPE(ref) == SVt_PVHV) {
inhv = (HV*)ref;
hv_iterinit(inhv);
HE *restrict e0 = hv_iternext(inhv);
if (!e0) { // empty hash: ambiguous shape -> empty result of the chosen/own shape
SV *restrict r = (want == FLT_AOH) ? newRV_noinc((SV*)newAV())
: newRV_noinc((SV*)newHV());
ST(0) = sv_2mortal(r);
XSRETURN(1);
}
NV O = obs_array[j];
SV**restrict col_key_sv = av_fetch(col_keys, j, 0);
hv_store_ent(expected_hv, *col_key_sv, newSVnv(E), 0);
stat += ((O - E) * (O - E)) / E;
}
expected_ref = newRV_noinc((SV*)expected_hv);
}
df = c - 1;
}
if (obs_matrix) {// Memory Cleanup for Matrices/Arrays
for (unsigned int i = 0; i < r; i++) {
safefree(obs_matrix[i]);
}
safefree(obs_matrix);
}
if (obs_array) safefree(obs_array);
if (row_keys) SvREFCNT_dec(row_keys);
if (col_keys) SvREFCNT_dec(col_keys);
NV p_val = get_p_value(stat, df);
// 3. Build the top-level results Hash (mimicking R's htest structure)
HV*restrict results = newHV();
HV*restrict statistic_hv = newHV();
hv_store(statistic_hv, "X-squared", 9, newSVnv(stat), 0);
hv_store(results, "statistic", 9, newRV_noinc((SV*)statistic_hv), 0);
HV*restrict parameter_hv = newHV();
hv_store(parameter_hv, "df", 2, newSViv(df), 0);
hv_store(results, "parameter", 9, newRV_noinc((SV*)parameter_hv), 0);
hv_store(results, "p.value", 7, newSVnv(p_val), 0);
hv_store(results, "expected", 8, expected_ref, 0);
hv_store(results, "observed", 8, SvREFCNT_inc(data_ref), 0);
if (input_type == SVt_PVAV) {
hv_store(results, "data.name", 9, newSVpv("Perl ArrayRef", 0), 0);
} else {
hv_store(results, "data.name", 9, newSVpv("Perl HashRef", 0), 0);
}
if (is_2d) {
if (yates) {
hv_store(results, "method", 6, newSVpv("Pearson's Chi-squared test with Yates' continuity correction", 0), 0);
} else {
hv_store(results, "method", 6, newSVpv("Pearson's Chi-squared test", 0), 0);
}
} else {
hv_store(results, "method", 6, newSVpv("Chi-squared test for given probabilities", 0), 0);
}
RETVAL = newRV_noinc((SV*)results);
}
OUTPUT:
RETVAL
PROTOTYPES: ENABLE
void write_table(...)
PPCODE:
{
SV *restrict data_sv = NULL;
SV *restrict file_sv = NULL;
unsigned int arg_idx = 0;
// Mimic the Perl shift logic
if (arg_idx < items && SvROK(ST(arg_idx))) {
int type = SvTYPE(SvRV(ST(arg_idx)));
if (type == SVt_PVHV || type == SVt_PVAV) {
data_sv = ST(arg_idx);
arg_idx++;
}
}
/* Only consume a positional file argument if it is a plain string that is
NOT one of the named option keys. Otherwise write_table(data=>..., file=>...)
would grab the literal string "data" as the filename.*/
if (arg_idx < items) {
SV *restrict cand = ST(arg_idx);
if (SvOK(cand) && !SvROK(cand)) {
const char *restrict k = SvPV_nolen(cand);
if (!(strEQ(k, "data") || strEQ(k, "file") || strEQ(k, "col.names") ||
strEQ(k, "row.names") || strEQ(k, "sep") || strEQ(k, "delim") ||
strEQ(k, "undef.val") || strEQ(k, "tex") ||
strEQ(k, "tex.col.align") || strEQ(k, "tex.size") ||
strEQ(k, "tex.comment") || strEQ(k, "tex.bold.1st.col") ||
strEQ(k, "tex.format") || strEQ(k, "tex.longtable") ||
strEQ(k, "tex.longtable.head") ||
strEQ(k, "xlsx") || strEQ(k, "xlsx.sheet") ||
strEQ(k, "xlsx.comment") || strEQ(k, "xlsx.freeze.rows") ||
strEQ(k, "xlsx.freeze.cols"))) {
file_sv = cand;
arg_idx++;
}
}
}
const char *restrict sep = ",";
bool explicit_sep = 0; // Track if delimiter was manually specified
/* default undef cells to a true empty value ("") instead of NULL.
With print_string_row emitting zero-length fields bare (no quotes), an
undef cell now prints as nothing at all: a,,c -- not a,'',c or a,"",c.
'undef.val' => 'NA' (etc.) still overrides this.*/
const char *restrict undef_val = "";
SV *restrict row_names_sv = sv_2mortal(newSViv(0));
SV *restrict col_names_sv = NULL;
/* LaTeX tabular output. 'tex' selects LaTeX for the main output file; the
remaining tex.* keys tune the rendering. tex_opt is tri-state: -1 = not
given (auto-detect from a ".tex" file name), 0 = off, 1 = on.*/
short int tex_opt = -1;
const char *restrict tex_align = "c"; // per-column alignment: c / l / r
const char *restrict tex_size = NULL; // optional size directive, e.g. \small
SV *restrict tex_comment = NULL; // string or array ref of % comment lines
bool tex_bold1 = 1; // bold the first column of each data row
bool tex_format = 0; // %.4g-format numeric cells
bool tex_longtable = 0; // body only, for \input into a longtable
/* Generate longtable's own repeat-header machinery (\endfirsthead / \endhead /
\endfoot) instead of a plain header row. A non-numeric value is the caption
used on continuation pages. Implies tex.longtable.*/
SV *restrict tex_longtable_head = NULL;
/* .xlsx (Excel) output, dependency-free. xlsx_opt is tri-state like tex_opt:
-1 = auto-detect from a ".xlsx" file name, 0 = off, 1 = on.*/
short int xlsx_opt = -1;
for (size_t i = 0; i < n; i++) {
if (rank_x[i] != floor(rank_x[i]) || rank_y[i] != floor(rank_y[i])) {
has_ties = 1;
break;
}
}
bool do_exact;
if (!exact_sv || !SvOK(exact_sv))
do_exact = (n < 10) && !has_ties;
else
do_exact = SvTRUE(exact_sv) ? 1 : 0;
if (do_exact) {
statistic = S_stat;
p_value = spearman_exact_pvalue(S_stat, n, alternative);
} else {
NV r = estimate;
/*NOTE: R silently ignores continuity correction for Spearman.
The adjustment below is non-standard; a warning is emitted
so callers are not silently misled.*/
if (continuity) {
warn("cor_test: continuity correction is not defined for Spearman in R and is ignored here");
}
NV denom_t = 1.0 - r * r;
if (denom_t <= 0.0)
statistic = (r > 0.0) ? INFINITY : -INFINITY;
else
statistic = r * sqrt((NV)(n - 2) / denom_t);
p_value = get_t_pvalue(statistic, (NV)(n - 2), alternative);
}
Safefree(rank_x); Safefree(rank_y);
} else {
Safefree(x); Safefree(y);
croak("Unknown method '%s': must be 'pearson', 'kendall', or 'spearman'", method);
}
Safefree(x); Safefree(y);
rhv = newHV();
hv_stores(rhv, "estimate", newSVnv(estimate));
hv_stores(rhv, "p.value", newSVnv(p_value));
hv_stores(rhv, "statistic", newSVnv(statistic));
hv_stores(rhv, "method", newSVpv(method, 0));
hv_stores(rhv, "alternative", newSVpv(alternative, 0));
if (is_pearson) {
hv_stores(rhv, "parameter", newSVnv(df));
AV *restrict ci_av = newAV();
av_push(ci_av, newSVnv(ci_lower));
av_push(ci_av, newSVnv(ci_upper));
hv_stores(rhv, "conf.int", newRV_noinc((SV*)ci_av));
}
RETVAL = newRV_noinc((SV*)rhv);
}
OUTPUT:
RETVAL
void shapiro_test(data)
SV *data
PREINIT:
AV *restrict av;
HV *restrict ret_hash;
size_t n_raw, n = 0;
NV *restrict x, w = 0.0, p_val = 0.0, mean = 0.0, ssq = 0.0;
PPCODE:
if (!SvROK(data) || SvTYPE(SvRV(data)) != SVt_PVAV) {
croak("Expected an array reference");
}
av = (AV *)SvRV(data);
n_raw = av_len(av) + 1;
Newx(x, n_raw, NV);
// Extract variables and calculate mean (skipping undefined/NaN values)
for (size_t i = 0; i < n_raw; i++) {
SV **restrict elem = av_fetch(av, i, 0);
if (elem && SvOK(*elem)) {
NV val = SvNV(*elem);
if (!isnan(val)) {
x[n] = val;
mean += val;
n++;
}
}
}
if (n < 3 || n > 5000) {
Safefree(x);
croak("Sample size must be between 3 and 5000 (R's limit)");
}
mean /= n;
for (size_t i = 0; i < n; i++) {// Calculate Sum of Squares
ssq += (x[i] - mean) * (x[i] - mean);
}
if (ssq == 0.0) {
Safefree(x);
croak("Data is perfectly constant; cannot compute Shapiro-Wilk test");
}
qsort(x, n, sizeof(NV), compare_doubles);
// --- Core AS R94 Algorithm: Weights and Statistic W
if (n == 3) {
NV a_val = 0.7071067811865475; // sqrt(1/2)
NV b_val = a_val * (x[2] - x[0]);
w = (b_val * b_val) / ssq;
if (w < 0.75) w = 0.75;
// Exact P-value for n=3
p_val = 1.90985931710274 * (asin(sqrt(w)) - 1.04719755119660);
} else {
NV *restrict m, *restrict a;
NV sum_m2 = 0.0, b_val = 0.0;
Newx(m, n, NV);
Newx(a, n, NV);
for (size_t i = 0; i < n; i++) {
m[i] = inverse_normal_cdf((i + 1.0 - 0.375) / (n + 0.25));
sum_m2 += m[i] * m[i];
}
NV u = 1.0 / sqrt((NV)n);
NV 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);
a[n-1] = a_n;
a[0] = -a_n;
if (n == 4 || n == 5) {
NV eps = (sum_m2 - 2.0 * m[n-1]*m[n-1]) / (1.0 - 2.0 * a_n*a_n);
for (unsigned int i = 1; i < n-1; i++) {
a[i] = m[i] / sqrt(eps);
}
} else {
NV a_n1 = -3.582633*pow(u,5) + 5.682633*pow(u,4) - 1.752461*pow(u,3) - 0.293762*pow(u,2) + 0.042981*u + m[n-2]/sqrt(sum_m2);
a[n-2] = a_n1;
q = (1.0 - gamma) * x[j] + gamma * x[j + 1];
}
// --- Format hash key with Epsilon guarding ---
char key[32];
double pct = (double)(p * 100.0); // Safe to cast to double just for formatting
double pct_rounded = floor(pct + 0.5); // C89 safe rounding
// Use 1e-9 epsilon check instead of strict integer equality
if (fabs(pct - pct_rounded) < 1e-9) {
snprintf(key, sizeof(key), "%.0f%%", pct_rounded);
} else {
snprintf(key, sizeof(key), "%.1f%%", pct);
}
hv_store(res_hv, key, strlen(key), newSVnv(q), 0);
}
Safefree(x); Safefree(probs);
RETVAL = newRV_noinc((SV*)res_hv);
}
OUTPUT:
RETVAL
NV mean(...)
PROTOTYPE: @
INIT:
NV total = 0;
size_t count = 0;
CODE:
for (size_t i = 0; i < items; i++) {
SV* restrict arg = ST(i);
if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
AV* restrict av = (AV*)SvRV(arg);
SSize_t len = av_len(av) + 1;
for (SSize_t j = 0; j < len; j++) {
SV** restrict tv = av_fetch(av, j, 0);
if (tv && SvOK(*tv)) {
total += SvNV(*tv);
count++;
} else {
croak("mean: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
} else if (SvOK(arg)) {
total += SvNV(arg);
count++;
} else {
croak("mean: undefined value at argument index %" UVuf, (UV)i);
}
}
if (count == 0) croak("mean needs >= 1 element");
RETVAL = total / count;
OUTPUT:
RETVAL
void mode(...)
PROTOTYPE: @
PREINIT:
HV *restrict counts;
HV *restrict originals;
size_t max_count = 0, arg_count = 0;
HE *restrict he;
PPCODE:
//counts: string(value) -> occurrence count
//originals: string(value) -> SV* first-seen original
counts = (HV *)sv_2mortal((SV *)newHV());
originals = (HV *)sv_2mortal((SV *)newHV());
for (size_t i = 0; i < items; i++) {
SV *restrict arg = ST(i);
if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
AV *restrict av = (AV *)SvRV(arg);
SSize_t len = av_len(av) + 1;
for (size_t j = 0; j < len; j++) {
SV **restrict tv = av_fetch(av, j, 0);
if (tv && SvOK(*tv)) {
STRLEN klen;
const char *restrict key = SvPV(*tv, klen);
SV **restrict slot = hv_fetch(counts, key, klen, 1);
if (!slot) croak("mode: internal hash error");
size_t cnt = SvOK(*slot) ? SvIV(*slot) + 1 : 1;
sv_setiv(*slot, cnt);
if (cnt > max_count) max_count = cnt;
if (cnt == 1)
hv_store(originals, key, klen, newSVsv(*tv), 0);
arg_count++;
} else {
croak("mode: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
} else if (SvOK(arg)) {
STRLEN klen;
const char *restrict key = SvPV(arg, klen);
SV **restrict slot = hv_fetch(counts, key, klen, 1);
if (!slot) croak("mode: internal hash error");
size_t cnt = SvOK(*slot) ? SvIV(*slot) + 1 : 1;
sv_setiv(*slot, cnt);
if (cnt > max_count) max_count = cnt;
if (cnt == 1)
hv_store(originals, key, klen, newSVsv(arg), 0);
arg_count++;
} else {
croak("mode: undefined value at argument index %" UVuf, (UV)i);
}
}
if (arg_count == 0)
croak("mode needs >= 1 element");
hv_iterinit(counts);
while ((he = hv_iternext(counts))) {
if (SvIV(hv_iterval(counts, he)) == max_count) {
STRLEN klen;
const char *restrict key = HePV(he, klen);
SV **restrict orig = hv_fetch(originals, key, klen, 0);
mXPUSHs(orig ? newSVsv(*orig) : newSVpvn(key, klen));
}
}
NV sum(...)
PROTOTYPE: @
INIT:
croak("sum: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
} else if (SvOK(arg)) {
total += SvNV(arg);
count++;
} else {
croak("sum: undefined value at argument index %" UVuf, (UV)i);
}
}
if (count == 0) croak("sum needs >= 1 element");
RETVAL = total;
OUTPUT:
RETVAL
NV sd(...)
PROTOTYPE: @
INIT:
NV mean = 0.0, M2 = 0.0;
size_t count = 0;
CODE:
for (size_t i = 0; i < items; i++) { // Single Pass Standard Deviation via Welford's Algorithm
SV* restrict arg = ST(i);
if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
AV* restrict av = (AV*)SvRV(arg);
SSize_t len = av_len(av) + 1;
for (size_t j = 0; j < len; j++) {
SV** restrict tv = av_fetch(av, j, 0);
if (tv && SvOK(*tv)) {
count++;
NV val = SvNV(*tv);
NV delta = val - mean;
mean += delta / count;
M2 += delta * (val - mean);
} else {
croak("sd: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
} else if (SvOK(arg)) {
count++;
NV val = SvNV(arg);
NV delta = val - mean;
mean += delta / count;
M2 += delta * (val - mean);
} else {
croak("sd: undefined value at argument index %" UVuf, (UV)i);
}
}
if (count < 2) croak("sd needs >= 2 elements");
RETVAL = sqrt(M2 / (count - 1));
OUTPUT:
RETVAL
void uniq(...)
PROTOTYPE: @
PREINIT:
HV*restrict seen;
AV*restrict out;
size_t n, k;
int gimme;
PPCODE:
n = 0;
gimme = GIMME_V;
seen = (HV*)sv_2mortal((SV*)newHV());
out = (AV*)sv_2mortal((SV*)newAV());
for (size_t i = 0; i < items; i++) {
SV* restrict arg = ST(i);
if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
AV* restrict av = (AV*)SvRV(arg);
size_t len = av_len(av) + 1;
for (size_t j = 0; j < len; j++) {
SV** restrict tv = av_fetch(av, j, 0);
if (tv && SvOK(*tv)) {
STRLEN klen;
const char*restrict key = SvPV(*tv, klen);
I32 hklen = SvUTF8(*tv) ? -(I32)klen : (I32)klen;
if (!hv_exists(seen, key, hklen)) {
(void)hv_store(seen, key, hklen, &PL_sv_undef, 0);
if (gimme != G_SCALAR)
av_push(out, newSVsv(*tv));
n++;
}
} else {
croak("uniq: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
} else if (SvOK(arg)) {
STRLEN klen;
const char*restrict key = SvPV(arg, klen);
I32 hklen = SvUTF8(arg) ? -(I32)klen : (I32)klen;
if (!hv_exists(seen, key, hklen)) {
(void)hv_store(seen, key, hklen, &PL_sv_undef, 0);
if (gimme != G_SCALAR)
av_push(out, newSVsv(arg));
n++;
}
} else {
croak("uniq: undefined value at argument index %" UVuf, (UV)i);
}
}
if (gimme == G_SCALAR) {
XPUSHs(sv_2mortal(newSVuv(n)));
} else {
size_t outlen = av_len(out) + 1;
EXTEND(SP, (SSize_t)outlen);
for (k = 0; k < outlen; k++)
PUSHs(sv_2mortal(av_shift(out)));
}
NV var(...)
PROTOTYPE: @
INIT:
NV mean = 0.0, M2 = 0.0;
size_t count = 0;
CODE:
// Single Pass Variance via Welford's Algorithm
for (size_t i = 0; i < items; i++) {
SV* restrict arg = ST(i);
if (SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV) {
AV* restrict av = (AV*)SvRV(arg);
size_t len = av_len(av) + 1;
(pow(stderr_x2, 2) / ((NV)nx - 1.0) + pow(stderr_y2, 2) / ((NV)ny - 1.0));
}
} else {
const size_t nx = t_test_scan(aTHX_ x_av, &mean_x, &var_x);
if (nx < 2) croak("t_test: 'x' needs at least 2 elements");
cint_est = mean_x;
std_err = sqrt(var_x / (NV)nx);
df = (NV)nx - 1.0;
constant_scale = fabs(mean_x);
}
/*R stops once the standard error has sunk into the rounding noise of
the data's own magnitude, not only when the variance is exactly zero.
An absolute test lets through a sample whose spread a double cannot
resolve at that scale and reports the noise as an enormous t: four
values around 1e10 differing by 1e-5 gave t = 4e15, p = 3e-47. The
exactly-zero case is R's NaN, croaked rather than returned.*/
if (std_err == 0.0
|| (isfinite(std_err) && std_err < 10.0 * DBL_EPSILON * constant_scale))
croak("t_test: data are essentially constant");
t_stat = (cint_est - mu) / std_err;
p_val = get_t_pvalue(t_stat, df, alternative);
HV*restrict results = newHV();
switch (estimates) {
case EST_MEAN_DIFF:
hv_store(results, "estimate", 8, newSVnv(cint_est), 0);
break;
case EST_BOTH:
hv_store(results, "estimate_x", 10, newSVnv(mean_x), 0);
hv_store(results, "estimate_y", 10, newSVnv(mean_y), 0);
break;
default:
hv_store(results, "estimate", 8, newSVnv(mean_x), 0);
}
NV alpha = 1.0 - conf_level, t_crit, ci_lower, ci_upper;
if (strcmp(alternative, "less") == 0) {
t_crit = qt_tail(df, alpha);
ci_lower = -INFINITY;
ci_upper = cint_est + t_crit * std_err;
} else if (strcmp(alternative, "greater") == 0) {
t_crit = qt_tail(df, alpha);
ci_lower = cint_est - t_crit * std_err;
ci_upper = INFINITY;
} else {
t_crit = qt_tail(df, alpha / 2.0);
ci_lower = cint_est - t_crit * std_err;
ci_upper = cint_est + t_crit * std_err;
}
AV*restrict conf_int = newAV();
av_push(conf_int, newSVnv(ci_lower));
av_push(conf_int, newSVnv(ci_upper));
hv_store(results, "statistic", 9, newSVnv(t_stat), 0);
hv_store(results, "df", 2, newSVnv(df), 0);
hv_store(results, "p_value", 7, newSVnv(p_val), 0);
hv_store(results, "conf_int", 8, newRV_noinc((SV*)conf_int), 0);
RETVAL = newRV_noinc((SV*)results);
}
OUTPUT:
RETVAL
void prop_test(...)
PPCODE:
{
/*Test of equality of proportions / a single proportion against a target.
Faithful port of R's stats::prop.test (Pearson chi-square on the 2xk
table of successes/failures, Yates correction for k<=2, Wilson score CI
for one proportion and a Wald CI for a difference of two).*/
if (items < 2)
croak("Usage: prop_test(\\@successes, \\@trials, p => ..., "
"alternative => 'two.sided', conf.level => 0.95, correct => 1)\n"
" or prop_test($x, $n, ...) for a single sample");
const char *restrict alt = "two.sided";
NV conf_level = 0.95;
bool correct = 1;
SV *restrict p_sv = NULL;
for (int i = 2; i + 1 < items; i += 2) {
const char *restrict key = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(key, "p")) p_sv = v;
else if (strEQ(key, "alternative")) alt = SvPV_nolen(v);
else if (strEQ(key, "conf_level") || strEQ(key, "conf.level")) conf_level = SvNV(v);
else if (strEQ(key, "correct")) correct = SvTRUE(v) ? 1 : 0;
else croak("prop_test: unknown argument '%s'", key);
}
if (!(conf_level > 0.0 && conf_level < 1.0))
croak("prop_test: conf.level must be between 0 and 1");
if (strNE(alt, "two.sided") && strNE(alt, "less") && strNE(alt, "greater"))
croak("prop_test: alternative must be 'two.sided', 'less' or 'greater'");
//--- read x (successes) and n (trials): each a scalar or an array ref ---
NV *restrict x = NULL, *restrict nn = NULL, *restrict pnull = NULL;
size_t k = 0;
{
SV *restrict xsv = ST(0), *restrict nsv = ST(1);
if (SvROK(xsv) && SvTYPE(SvRV(xsv)) == SVt_PVAV) {
AV *av = (AV*)SvRV(xsv); k = (size_t)(av_len(av) + 1);
if (k == 0) croak("prop_test: 'x' is empty");
Newx(x, k, NV);
for (size_t i = 0; i < k; i++) { SV **e = av_fetch(av, i, 0); x[i] = (e && *e) ? SvNV(*e) : NAN; }
} else { k = 1; Newx(x, 1, NV); x[0] = SvNV(xsv); }
size_t kn;
if (SvROK(nsv) && SvTYPE(SvRV(nsv)) == SVt_PVAV) {
AV *restrict av = (AV*)SvRV(nsv); kn = (size_t)(av_len(av) + 1);
Newx(nn, kn, NV);
for (size_t i = 0; i < kn; i++) { SV **e = av_fetch(av, i, 0); nn[i] = (e && *e) ? SvNV(*e) : NAN; }
} else { kn = 1; Newx(nn, 1, NV); nn[0] = SvNV(nsv); }
if (kn != k) { Safefree(x); Safefree(nn); croak("prop_test: 'x' and 'n' must have the same length"); }
}
for (size_t i = 0; i < k; i++) {
if (nn[i] <= 0) { Safefree(x); Safefree(nn); croak("prop_test: elements of 'n' must be positive"); }
if (x[i] < 0) { Safefree(x); Safefree(nn); croak("prop_test: elements of 'x' must be nonnegative"); }
if (x[i] > nn[i]) { Safefree(x); Safefree(nn); croak("prop_test: elements of 'x' must not exceed 'n'"); }
}
//--- null probabilities and degrees of freedom ---
bool p_is_null; //true => testing equality (pooled p, df = k-1)
Newx(pnull, k, NV);
if (p_sv && SvOK(p_sv)) {
p_is_null = FALSE;
if (SvROK(p_sv) && SvTYPE(SvRV(p_sv)) == SVt_PVAV) {
AV *restrict av = (AV*)SvRV(p_sv);
if ((size_t)(av_len(av) + 1) != k) { Safefree(x); Safefree(nn); Safefree(pnull); croak("prop_test: 'p' must have the same length as 'x'"); }
for (size_t i = 0; i < k; i++) { SV **e = av_fetch(av, i, 0); pnull[i] = (e && *e) ? SvNV(*e) : NAN; }
else { ci_lo = -1.0; ci_hi = (delta + width > 1.0) ? 1.0 : delta + width; }
have_ci = TRUE;
}
int df = p_is_null ? (int)(k - 1) : (int)k;
//--- Pearson chi-square with (capped) Yates correction ---
NV stat = 0.0;
for (size_t i = 0; i < k; i++) {
NV E0 = nn[i] * pnull[i], E1 = nn[i] * (1.0 - pnull[i]);
NV o0 = x[i], o1 = nn[i] - x[i];
NV d0 = fabs(o0 - E0) - YATES; //R does not floor this at 0
NV d1 = fabs(o1 - E1) - YATES;
stat += d0 * d0 / E0 + d1 * d1 / E1;
}
NV p_value;
if (strEQ(alt, "two.sided")) {
p_value = get_p_value(stat, df);
} else {
NV z = (k == 1) ? ((est[0] > pnull[0]) - (est[0] < pnull[0])) * sqrt(stat)
: ((delta > 0.0) - (delta < 0.0)) * sqrt(stat);
p_value = strEQ(alt, "less") ? approx_pnorm(z) : 1.0 - approx_pnorm(z);
}
// Chi-square approximation warning, as in R
for (size_t i = 0; i < k; i++)
if (nn[i] * pnull[i] < 5.0 || nn[i] * (1.0 - pnull[i]) < 5.0) {
warn("prop_test: Chi-squared approximation may be incorrect");
break;
}
HV *restrict ret = newHV();
hv_stores(ret, "statistic", newSVnv(stat));
hv_stores(ret, "parameter", newSViv(df));
hv_stores(ret, "p_value", newSVnv(p_value));
hv_stores(ret, "alternative", newSVpv(alt, 0));
hv_stores(ret, "conf_level", newSVnv(conf_level));
{
char method[96];
if (k == 1) snprintf(method, sizeof method, "1-sample proportions test %s continuity correction", YATES > 0.0 ? "with" : "without");
else snprintf(method, sizeof method, "%zu-sample test for %s proportions %s continuity correction",
k, p_is_null ? "equality of" : "given", YATES > 0.0 ? "with" : "without");
hv_stores(ret, "method", newSVpv(method, 0));
}
{
AV *restrict ev = newAV();
for (size_t i = 0; i < k; i++) av_push(ev, newSVnv(est[i]));
hv_stores(ret, "estimate", newRV_noinc((SV*)ev));
}
if (have_ci) {
AV *restrict ci = newAV(); av_push(ci, newSVnv(ci_lo)); av_push(ci, newSVnv(ci_hi));
hv_stores(ret, "conf.int", newRV_noinc((SV*)ci));
}
Safefree(x); Safefree(nn); Safefree(pnull); Safefree(est);
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void mcnemar_test(...)
PPCODE:
{
/*McNemar's test for paired categorical data. Faithful port of R's
stats::mcnemar.test (chi-square on the off-diagonal disagreement,
with Yates continuity correction for a 2x2 table). An `exact => 1`
option gives the two-sided exact binomial test for a 2x2 table.*/
if (items < 1)
croak("Usage: mcnemar_test([[a,b],[c,d]], correct => 1, exact => 0)\n"
" or mcnemar_test(\\@x, \\@y, ...) # paired observations");
int correct = 1, exact = 0, opt_start;
size_t r = 0;
NV *restrict tab = NULL;
bool is_matrix = FALSE;
if (SvROK(ST(0)) && SvTYPE(SvRV(ST(0))) == SVt_PVAV
&& av_len((AV*)SvRV(ST(0))) >= 0) {
SV **restrict e0 = av_fetch((AV*)SvRV(ST(0)), 0, 0);
is_matrix = e0 && *e0 && SvROK(*e0) && SvTYPE(SvRV(*e0)) == SVt_PVAV;
}
if (is_matrix) {
AV *restrict m = (AV*)SvRV(ST(0));
r = (size_t)(av_len(m) + 1);
if (r < 2) croak("mcnemar_test: matrix must have at least two rows");
Newxz(tab, r * r, NV);
for (size_t i = 0; i < r; i++) {
SV **restrict row = av_fetch(m, i, 0);
if (!row || !*row || !SvROK(*row) || SvTYPE(SvRV(*row)) != SVt_PVAV)
{ Safefree(tab); croak("mcnemar_test: row %" UVuf " is not an array ref", (UV)i); }
AV *rv = (AV*)SvRV(*row);
if ((size_t)(av_len(rv) + 1) != r) { Safefree(tab); croak("mcnemar_test: matrix must be square"); }
for (size_t j = 0; j < r; j++) {
SV **restrict c = av_fetch(rv, j, 0);
NV v = (c && *c) ? SvNV(*c) : 0.0;
if (v < 0 || isnan(v)) { Safefree(tab); croak("mcnemar_test: entries must be nonnegative and finite"); }
tab[i * r + j] = v;
}
}
opt_start = 1;
} else {
if (items < 2 || !SvROK(ST(0)) || !SvROK(ST(1))
|| SvTYPE(SvRV(ST(0))) != SVt_PVAV || SvTYPE(SvRV(ST(1))) != SVt_PVAV)
croak("mcnemar_test: expected a square matrix or two array refs");
AV *restrict xa = (AV*)SvRV(ST(0)), *ya = (AV*)SvRV(ST(1));
size_t n = (size_t)(av_len(xa) + 1);
if ((size_t)(av_len(ya) + 1) != n) croak("mcnemar_test: 'x' and 'y' must have the same length");
//collect sorted unique levels across both vectors
char **restrict lev = NULL; size_t nlev = 0, cap = 8; Newx(lev, cap, char*);
for (size_t src = 0; src < 2; src++) {
AV *a = src ? ya : xa;
for (size_t i = 0; i < n; i++) {
SV **e = av_fetch(a, i, 0);
if (!e || !*e || !SvOK(*e)) continue;
STRLEN l; const char *s = SvPV(*e, l);
bool found = FALSE;
for (size_t k = 0; k < nlev; k++) if (strEQ(lev[k], s)) { found = TRUE; break; }
if (!found) { if (nlev >= cap) { cap *= 2; Renew(lev, cap, char*); } lev[nlev++] = savepvn(s, l); }
}
}
//sort levels lexically for a deterministic table order
for (size_t a = 0; a + 1 < nlev; a++) for (size_t b = a + 1; b < nlev; b++)
for (size_t k = 0; k < r; k++) { if (strEQ(lev[k], sx)) ix = k; if (strEQ(lev[k], sy)) iy = k; }
tab[ix * r + iy] += 1.0;
}
for (size_t k = 0; k < nlev; k++) Safefree(lev[k]);
Safefree(lev);
opt_start = 2;
}
for (int i = opt_start; i + 1 < items; i += 2) {
const char *restrict k = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(k, "correct")) correct = SvTRUE(v) ? 1 : 0;
else if (strEQ(k, "exact")) exact = SvTRUE(v) ? 1 : 0;
else { Safefree(tab); croak("mcnemar_test: unknown argument '%s'", k); }
}
if (exact && r != 2) { Safefree(tab); croak("mcnemar_test: exact test requires a 2x2 table"); }
HV *restrict ret = newHV();
if (exact) {// two-sided exact binomial test of the discordant pairs, b ~ Bin(b+c, 0.5)
NV b = tab[0 * 2 + 1], c = tab[1 * 2 + 0];
NV nn = b + c, m = (b < c) ? b : c;
NV p_value;
if (nn == 0.0) p_value = 1.0;
else {
NV s = 0.0, lognn2 = nn * log(2.0);
for (NV kk = 0.0; kk <= m; kk += 1.0) s += exp(ft_lchoose((long)nn, (long)kk) - lognn2);
p_value = 2.0 * s; if (p_value > 1.0) p_value = 1.0;
}
hv_stores(ret, "statistic", newSVnv(b));
hv_stores(ret, "p_value", newSVnv(p_value));
hv_stores(ret, "method", newSVpv("McNemar's test (exact binomial)", 0));
} else {
bool use_cc = 0;
if (correct && r == 2) {
for (size_t i = 0; i < r && !use_cc; i++)
for (size_t j = 0; j < r; j++)
if (tab[i * r + j] != tab[j * r + i]) { use_cc = 1; break; }
}
NV stat = 0.0;
for (size_t i = 0; i < r; i++)
for (size_t j = i + 1; j < r; j++) {
NV diff = tab[i * r + j] - tab[j * r + i];
NV sum = tab[i * r + j] + tab[j * r + i];
if (sum <= 0.0) continue;
NV num = use_cc ? (fabs(diff) - 1.0) : diff;
stat += num * num / sum;
}
int df = (int)(r * (r - 1) / 2);
NV p_value = get_p_value(stat, df);
hv_stores(ret, "statistic", newSVnv(stat));
hv_stores(ret, "parameter", newSViv(df));
hv_stores(ret, "p_value", newSVnv(p_value));
hv_stores(ret, "method", newSVpv(use_cc ?
"McNemar's Chi-squared test with continuity correction" :
"McNemar's Chi-squared test", 0));
}
Safefree(tab);
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void dunn_test(...)
PPCODE:
{
/*Dunn's (1964) post-hoc test following a Kruskal-Wallis test: pairwise
rank-mean comparisons using the shared ranking and tie correction, with
a family-wise / FDR adjustment. Two-sided p-values (as in FSA::dunnTest).
Validated against the canonical formula implemented in base R.*/
if (items < 2 || !SvROK(ST(0)) || !SvROK(ST(1))
|| SvTYPE(SvRV(ST(0))) != SVt_PVAV || SvTYPE(SvRV(ST(1))) != SVt_PVAV)
croak("Usage: dunn_test(\\@values, \\@groups, method => 'holm')");
const char *restrict method = "holm";
for (int i = 2; i + 1 < items; i += 2) {
const char *key = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(key, "method")) method = SvPV_nolen(v);
else croak("dunn_test: unknown argument '%s'", key);
}
char meth[32]; strncpy(meth, method, 31); meth[31] = '\0';
for (unsigned i = 0; meth[i]; i++) meth[i] = tolower(meth[i]);
if (strEQ(meth, "fdr")) strcpy(meth, "bh");
if (strEQ(meth, "holm-sidak")) strcpy(meth, "hs");
AV *restrict xa = (AV*)SvRV(ST(0)), *ga = (AV*)SvRV(ST(1));
size_t raw = (size_t)(av_len(xa) + 1);
if ((size_t)(av_len(ga) + 1) != raw) croak("dunn_test: values and groups must have the same length");
//gather complete (value, group) pairs
NV *restrict x = NULL; char **restrict glab = NULL;
Newx(x, raw, NV); Newx(glab, raw, char*);
size_t N = 0;
for (size_t i = 0; i < raw; i++) {
SV **xv = av_fetch(xa, i, 0), **gv = av_fetch(ga, i, 0);
if (!xv || !*xv || !SvOK(*xv) || !looks_like_number(*xv)) continue;
if (!gv || !*gv || !SvOK(*gv)) continue;
NV val = SvNV(*xv); if (isnan(val)) continue;
STRLEN l; const char *s = SvPV(*gv, l);
x[N] = val; glab[N] = savepvn(s, l); N++;
}
if (N < 3) { for (size_t i = 0; i < N; i++) Safefree(glab[i]); Safefree(x); Safefree(glab); croak("dunn_test: not enough complete observations"); }
// sorted unique group levels
char **restrict lev = NULL; size_t k = 0, cap = 8; Newx(lev, cap, char*);
for (size_t i = 0; i < N; i++) {
bool found = FALSE;
for (size_t j = 0; j < k; j++) if (strEQ(lev[j], glab[i])) { found = TRUE; break; }
if (!found) { if (k >= cap) { cap *= 2; Renew(lev, cap, char*); } lev[k++] = savepv(glab[i]); }
}
for (size_t a = 0; a + 1 < k; a++) for (size_t b = a + 1; b < k; b++)
if (strcmp(lev[a], lev[b]) > 0) { char *t = lev[a]; lev[a] = lev[b]; lev[b] = t; }
if (k < 2) { for (size_t i = 0; i < N; i++) Safefree(glab[i]); for (size_t j = 0; j < k; j++) Safefree(lev[j]); Safefree(x); Safefree(glab); Safefree(lev); croak("dunn_test: need at least two groups"); }
//ranks over all observations (tie-averaged)
NV *restrict r = NULL; Newx(r, N, NV);
rank_data(x, r, N);
//per-group rank sums and sizes
NV *restrict rsum = NULL; size_t *restrict ns = NULL;
Newxz(rsum, k, NV); Newxz(ns, k, size_t);
for (size_t i = 0; i < N; i++) {
size_t gi = 0;
for (size_t j = 0; j < k; j++) if (strEQ(lev[j], glab[i])) { gi = j; break; }
rsum[gi] += r[i]; ns[gi]++;
}
//tie correction: sum over distinct values of (t^3 - t)
NV *restrict xs = NULL; Newx(xs, N, NV);
memcpy(xs, x, N * sizeof(NV));
qsort(xs, N, sizeof(NV), cmp_nv3);
NV tsum = 0.0;
{
size_t a = 0;
while (a < N) {
size_t b = a;
while (b + 1 < N && xs[b + 1] == xs[a]) b++;
NV t = (NV)(b - a + 1);
tsum += t * t * t - t;
a = b + 1;
}
}
Safefree(xs);
NV Nf = (NV)N;
NV sigma_base = (Nf * (Nf + 1.0)) / 12.0 - tsum / (12.0 * (Nf - 1.0));
size_t m = k * (k - 1) / 2;
NV *restrict z = NULL, *restrict praw = NULL, *restrict padj = NULL;
Newx(z, m, NV); Newx(praw, m, NV); Newx(padj, m, NV);
size_t *restrict gi_ = NULL, *restrict gj_ = NULL;
Newx(gi_, m, size_t); Newx(gj_, m, size_t);
size_t c = 0;
for (size_t i = 0; i < k; i++)
for (size_t j = i + 1; j < k; j++) {
NV rbar_i = rsum[i] / ns[i], rbar_j = rsum[j] / ns[j];
NV se = sqrt(sigma_base * (1.0 / ns[i] + 1.0 / ns[j]));
NV zz = (rbar_i - rbar_j) / se;
z[c] = zz;
praw[c] = 2.0 * (1.0 - approx_pnorm(fabs(zz)));
if (praw[c] > 1.0) praw[c] = 1.0;
gi_[c] = i; gj_[c] = j;
c++;
}
dunn_padjust(praw, m, meth, padj);
AV *restrict out = newAV();
for (size_t t = 0; t < m; t++) {
HV *restrict h = newHV();
char comp[256];
snprintf(comp, sizeof comp, "%s - %s", lev[gi_[t]], lev[gj_[t]]);
hv_stores(h, "comparison", newSVpv(comp, 0));
hv_stores(h, "group1", newSVpv(lev[gi_[t]], 0));
hv_stores(h, "group2", newSVpv(lev[gj_[t]], 0));
hv_stores(h, "Z", newSVnv(z[t]));
hv_stores(h, "p_value", newSVnv(praw[t]));
hv_stores(h, "p_adjust", newSVnv(padj[t]));
av_push(out, newRV_noinc((SV*)h));
}
for (size_t i = 0; i < N; i++) Safefree(glab[i]);
for (size_t j = 0; j < k; j++) Safefree(lev[j]);
Safefree(x); Safefree(glab); Safefree(lev); Safefree(r);
Safefree(rsum); Safefree(ns); Safefree(z); Safefree(praw); Safefree(padj);
Safefree(gi_); Safefree(gj_);
ST(0) = sv_2mortal(newRV_noinc((SV*)out));
XSRETURN(1);
}
void friedman_test(...)
PPCODE:
{
/*Friedman rank-sum test for an unreplicated complete block design.
Input is a matrix (array of array refs) with one block/subject per row
and one treatment/condition per column. Faithful port of R's
stats::friedman.test, including the tie correction.*/
if (items < 1 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV)
croak("Usage: friedman_test([[..row1..],[..row2..], ...]) # rows = blocks, cols = treatments");
AV *restrict m = (AV*)SvRV(ST(0));
size_t nrow_raw = (size_t)(av_len(m) + 1);
if (nrow_raw < 2) croak("friedman_test: need at least two blocks (rows)");
// determine k from the first row
SV **restrict r0 = av_fetch(m, 0, 0);
if (!r0 || !*r0 || !SvROK(*r0) || SvTYPE(SvRV(*r0)) != SVt_PVAV)
croak("friedman_test: each row must be an array ref");
size_t k = (size_t)(av_len((AV*)SvRV(*r0)) + 1);
if (k < 2) croak("friedman_test: need at least two treatments (columns)");
NV *restrict colsum = NULL; Newxz(colsum, k, NV);
NV *restrict rowbuf = NULL; Newx(rowbuf, k, NV);
NV *restrict ranks = NULL; Newx(ranks, k, NV);
NV *restrict sorted = NULL; Newx(sorted, k, NV);
NV tie_sum = 0.0;
size_t n = 0; // complete blocks actually used
for (size_t i = 0; i < nrow_raw; i++) {
SV **restrict rr = av_fetch(m, i, 0);
if (!rr || !*rr || !SvROK(*rr) || SvTYPE(SvRV(*rr)) != SVt_PVAV)
{ Safefree(colsum); Safefree(rowbuf); Safefree(ranks); Safefree(sorted); croak("friedman_test: row %" UVuf " is not an array ref", (UV)i); }
AV *rv = (AV*)SvRV(*rr);
if ((size_t)(av_len(rv) + 1) != k)
{ Safefree(colsum); Safefree(rowbuf); Safefree(ranks); Safefree(sorted); croak("friedman_test: all rows must have the same number of columns"); }
bool complete = TRUE;
for (size_t j = 0; j < k; j++) {
SV **c = av_fetch(rv, j, 0);
if (!c || !*c || !SvOK(*c) || !looks_like_number(*c)) { complete = FALSE; break; }
rowbuf[j] = SvNV(*c);
if (isnan(rowbuf[j])) { complete = FALSE; break; }
}
if (!complete) continue; //drop incomplete blocks, like R's complete.cases
rank_data(rowbuf, ranks, k);
for (size_t j = 0; j < k; j++) colsum[j] += ranks[j];
//tie correction: sum over tie groups of (u^3 - u) within this block
memcpy(sorted, rowbuf, k * sizeof(NV));
qsort(sorted, k, sizeof(NV), cmp_nv3);
size_t a = 0;
while (a < k) {
size_t b = a;
while (b + 1 < k && sorted[b + 1] == sorted[a]) b++;
NV u = (NV)(b - a + 1);
tie_sum += u * u * u - u;
a = b + 1;
}
n++;
}
Safefree(rowbuf); Safefree(ranks); Safefree(sorted);
if (n < 1) { Safefree(colsum); croak("friedman_test: no complete blocks"); }
NV nf = (NV)n, kf = (NV)k;
NV ssq = 0.0, target = nf * (kf + 1.0) / 2.0;
for (size_t j = 0; j < k; j++) { NV d = colsum[j] - target; ssq += d * d; }
Safefree(colsum);
NV denom = nf * kf * (kf + 1.0) - tie_sum / (kf - 1.0);
NV stat = 12.0 * ssq / denom;
int df = (int)(k - 1);
NV p_value = get_p_value(stat, df);
HV *ret = newHV();
hv_stores(ret, "statistic", newSVnv(stat));
hv_stores(ret, "parameter", newSViv(df));
hv_stores(ret, "p_value", newSVnv(p_value));
hv_stores(ret, "n", newSViv((int)n));
hv_stores(ret, "method", newSVpv("Friedman rank sum test", 0));
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void epi_2x2(...)
PPCODE:
{
NV a, b, c, d, conf_level = 0.95;
int correct = 0, opt_start;
if (items < 1)
croak("Usage: epi_2x2(a, b, c, d, conf_level => 0.95, correct => 0)\n"
" or epi_2x2([[a,b],[c,d]], ...) # rows=exposure, cols=outcome");
if (SvROK(ST(0))) {
epi_read_2x2(aTHX_ ST(0), "epi_2x2", &a, &b, &c, &d);
opt_start = 1;
} else {
if (items < 4)
croak("epi_2x2: need 4 cell counts a,b,c,d (or a single 2x2 array ref)");
a = SvNV(ST(0)); b = SvNV(ST(1)); c = SvNV(ST(2)); d = SvNV(ST(3));
if (a < 0 || b < 0 || c < 0 || d < 0)
croak("epi_2x2: cell counts must be non-negative");
opt_start = 4;
}
for (int i = opt_start; i + 1 < items; i += 2) {
const char *k = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(k, "conf_level") || strEQ(k, "conf.level")) conf_level = SvNV(v);
else if (strEQ(k, "correct")) correct = SvTRUE(v) ? 1 : 0;
else croak("epi_2x2: unknown argument '%s'", k);
}
if (!(conf_level > 0.0 && conf_level < 1.0))
croak("epi_2x2: conf_level must be between 0 and 1");
/*Haldane-Anscombe +0.5 when asked, or forced by a zero cell (avoids
division by zero / log of zero). Point estimates then shift too, so
the applied flag is reported back.*/
NV A = a, B = b, C = c, D = d; int corrected = 0;
if (correct || a == 0 || b == 0 || c == 0 || d == 0) {
A += 0.5; B += 0.5; C += 0.5; D += 0.5; corrected = 1;
}
NV z = inverse_normal_cdf(1.0 - (1.0 - conf_level) / 2.0);
NV n1 = A + B, n0 = C + D;
NV or_ = (A * D) / (B * C);
NV se_lor = sqrt(1.0 / A + 1.0 / B + 1.0 / C + 1.0 / D);
NV or_lo = or_ * exp(-z * se_lor), or_hi = or_ * exp(z * se_lor);
NV p1 = A / n1, p0 = C / n0;
NV rr = p1 / p0; //Katz log-RR variance
NV se_lrr = sqrt(B / (A * n1) + D / (C * n0));
NV rr_lo = rr * exp(-z * se_lrr), rr_hi = rr * exp(z * se_lrr);
NV rd = p1 - p0;
NV se_rd = sqrt(p1 * (1.0 - p1) / n1 + p0 * (1.0 - p0) / n0);
NV rd_lo = rd - z * se_rd, rd_hi = rd + z * se_rd;
HV *restrict ret = newHV();
hv_stores(ret, "method", newSVpv("2x2 epidemiological measures (Wald)", 0));
hv_stores(ret, "conf_level", newSVnv(conf_level));
hv_stores(ret, "correction", newSViv(corrected));
hv_stores(ret, "odds_ratio", newSVnv(or_));
hv_stores(ret, "risk_ratio", newSVnv(rr));
hv_stores(ret, "risk_diff", newSVnv(rd));
hv_stores(ret, "risk_exposed", newSVnv(p1));
hv_stores(ret, "risk_unexposed", newSVnv(p0));
hv_stores(ret, "nnt", newSVnv(1.0 / fabs(rd)));
#define EPI_CI(name, lo, hi) do { AV *ci = newAV(); \
av_push(ci, newSVnv(lo)); av_push(ci, newSVnv(hi)); \
hv_stores(ret, name, newRV_noinc((SV *)ci)); } while (0)
EPI_CI("odds_ratio_ci", or_lo, or_hi);
EPI_CI("risk_ratio_ci", rr_lo, rr_hi);
EPI_CI("risk_diff_ci", rd_lo, rd_hi);
#undef EPI_CI
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void cmh_test(...)
PPCODE:
{
if (items < 1 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV)
croak("Usage: cmh_test([ [a,b,c,d], [a,b,c,d], ... ], "
"conf_level => 0.95, correct => 1)");
AV *restrict strata = (AV *)SvRV(ST(0));
SSize_t K = av_len(strata) + 1;
if (K < 1) croak("cmh_test: need at least one 2x2 stratum");
NV conf_level = 0.95; int correct = 1;
for (int i = 1; i + 1 < items; i += 2) {
const char *k = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(k, "conf_level") || strEQ(k, "conf.level")) conf_level = SvNV(v);
else if (strEQ(k, "correct")) correct = SvTRUE(v) ? 1 : 0;
else croak("cmh_test: unknown argument '%s'", k);
}
if (!(conf_level > 0.0 && conf_level < 1.0))
croak("cmh_test: conf_level must be between 0 and 1");
NV sum_a = 0, E = 0, V = 0; //CMH statistic pieces
NV sumR = 0, sumS = 0; //Mantel-Haenszel common-OR num/den
NV vR = 0, vRS = 0, vS = 0; //Robins-Breslow-Greenland variance
for (SSize_t s = 0; s < K; s++) {
SV **restrict ep = av_fetch(strata, s, 0);
NV a, b, c, d;
epi_read_2x2(aTHX_ (ep ? *ep : &PL_sv_undef), "cmh_test", &a, &b, &c, &d);
NV n = a + b + c + d;
if (n <= 0) continue;
NV r1 = a + b, r2 = c + d, c1 = a + c, c2 = b + d;
sum_a += a;
E += r1 * c1 / n;
if (n > 1.0) V += (r1 * r2 * c1 * c2) / (n * n * (n - 1.0));
NV Rk = a * d / n, Sk = b * c / n;
sumR += Rk; sumS += Sk;
NV Pk = (a + d) / n, Qk = (b + c) / n;
vR += Pk * Rk;
vRS += Pk * Sk + Qk * Rk;
vS += Qk * Sk;
}
NV diff = fabs(sum_a - E) - (correct ? 0.5 : 0.0);
if (diff < 0) diff = 0;
NV chi = (V > 0) ? (diff * diff) / V : 0.0;
NV pval = get_p_value(chi, 1);
NV or_mh = (sumS > 0) ? sumR / sumS : NAN;
NV var_lnor = vR / (2.0 * sumR * sumR)
+ vRS / (2.0 * sumR * sumS)
+ vS / (2.0 * sumS * sumS);
NV z = inverse_normal_cdf(1.0 - (1.0 - conf_level) / 2.0);
NV or_lo = or_mh * exp(-z * sqrt(var_lnor)), or_hi = or_mh * exp(z * sqrt(var_lnor));
HV *ret = newHV();
hv_stores(ret, "method", newSVpv(correct
? "Mantel-Haenszel chi-squared test with continuity correction"
: "Mantel-Haenszel chi-squared test", 0));
hv_stores(ret, "statistic", newSVnv(chi));
hv_stores(ret, "parameter", newSViv(1)); //degrees of freedom
hv_stores(ret, "p_value", newSVnv(pval));
hv_stores(ret, "estimate", newSVnv(or_mh)); //common odds ratio
hv_stores(ret, "conf_level", newSVnv(conf_level));
hv_stores(ret, "correction", newSViv(correct));
hv_stores(ret, "k", newSViv((IV)K));
if (have_cutoff || have_frac) {
/*Binarize a continuous truth column, then DeLong on the pos/neg
score vectors -- same shape as bedroc's cutoff/active_frac.*/
SSize_t Ns = av_len(sav) + 1;
if (Ns != av_len(lav) + 1)
croak("auroc: y_true and y_score must be the same length");
if (Ns < 1) croak("auroc: need at least one observation");
size_t N = (size_t)Ns;
int *restrict act; Newxz(act, N, int);
if (have_frac) {
size_t n_a = (size_t)ceil(active_frac * (NV)N);
if (n_a < 1) n_a = 1;
if (N >= 2 && n_a > N - 1) n_a = N - 1;
NVIdx *restrict tmp; Newx(tmp, N, NVIdx);
for (size_t i = 0; i < N; i++) {
SV **restrict lp = av_fetch(lav, i, 0);
tmp[i].v = (lp && *lp) ? SvNV(*lp) : NAN; tmp[i].i = i;
}
qsort(tmp, N, sizeof(NVIdx), nvidx_cmp_asc);
if (frac_low) for (size_t k = 0; k < n_a; k++) act[tmp[k].i] = 1;
else for (size_t k = N - n_a; k < N; k++) act[tmp[k].i] = 1;
Safefree(tmp);
} else {
for (size_t i = 0; i < N; i++) {
SV **restrict lp = av_fetch(lav, i, 0);
act[i] = (((lp && *lp) ? SvNV(*lp) : NAN) >= cutoff);
}
}
NV *restrict P; Newx(P, N, NV); NV *restrict Q; Newx(Q, N, NV);
size_t m = 0, n = 0;
for (size_t i = 0; i < N; i++) {
SV **restrict sp = av_fetch(sav, i, 0);
NV s = (sp && *sp) ? SvNV(*sp) : NAN;
if (lower_pos) s = -s;
if (act[i]) P[m++] = s; else Q[n++] = s;
}
Safefree(act);
if (m == 0 || n == 0) {
Safefree(P); Safefree(Q);
croak("auroc: need both positive and negative labels%s",
have_cutoff ? " (check cutoff)" : "");
}
NV a, se; roc_delong(aTHX_ P, m, Q, n, &a, &se);
Safefree(P); Safefree(Q);
RETVAL = a;
} else {
/*Pre-built 0/1 (or string) labels: roc_split does the parallel-array
walk; note swapped args so labels are ST(0), scores ST(1).*/
NV *pos, *neg; size_t m, n;
roc_split(aTHX_ sav, lav, positive, lower_pos, &pos, &m, &neg, &n, "auroc");
NV a, se; roc_delong(aTHX_ pos, m, neg, n, &a, &se);
Safefree(pos); Safefree(neg);
RETVAL = a;
}
}
OUTPUT:
RETVAL
void roc(...)
PPCODE:
{
if (items < 2 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV
|| !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV)
croak("Usage: roc(\\@scores, \\@labels, positive => 1, "
"conf_level => 0.95, direction => '>')");
const char *restrict positive = "1"; NV conf_level = 0.95; int lower_pos = 0;
for (int i = 2; i + 1 < items; i += 2) {
const char *restrict k = SvPV_nolen(ST(i)); SV *restrict v = ST(i + 1);
if (strEQ(k, "positive")) positive = SvPV_nolen(v);
else if (strEQ(k, "conf_level") || strEQ(k, "conf.level")) conf_level = SvNV(v);
else if (strEQ(k, "direction")) { const char *restrict d = SvPV_nolen(v); lower_pos = (d[0] == '<'); }
else croak("roc: unknown argument '%s'", k);
}
if (!(conf_level > 0.0 && conf_level < 1.0))
croak("roc: conf_level must be between 0 and 1");
NV *pos, *neg; size_t m, n;
roc_split(aTHX_ (AV *)SvRV(ST(0)), (AV *)SvRV(ST(1)), positive, lower_pos,
&pos, &m, &neg, &n, "roc");
NV auc_val, se; roc_delong(aTHX_ pos, m, neg, n, &auc_val, &se);
NV z = inverse_normal_cdf(1.0 - (1.0 - conf_level) / 2.0);
NV lo = auc_val - z * se, hi = auc_val + z * se;
if (lo < 0.0) lo = 0.0; if (hi > 1.0) hi = 1.0;
// Curve + Youden by sweeping thresholds high -> low over all points
size_t N = m + n;
ROCPt *restrict pts; Newx(pts, N, ROCPt);
for (size_t i = 0; i < m; i++) { pts[i].score = pos[i]; pts[i].lab = 1; }
for (size_t j = 0; j < n; j++) { pts[m + j].score = neg[j]; pts[m + j].lab = 0; }
Safefree(pos); Safefree(neg);
qsort(pts, N, sizeof(ROCPt), rocpt_cmp_desc);
AV *restrict curve = newAV();
{ //leading operating point: threshold = +inf, nothing called positive
HV *restrict p0 = newHV();
hv_stores(p0, "threshold", newSVnv(INFINITY));
hv_stores(p0, "sensitivity", newSVnv(0.0));
hv_stores(p0, "specificity", newSVnv(1.0));
av_push(curve, newRV_noinc((SV *)p0));
}
NV TP = 0.0, FP = 0.0, P = (NV)m, Nn = (NV)n;
NV best_j = -2.0, best_thr = 0.0, best_sens = 0.0, best_spec = 0.0;
size_t i = 0;
while (i < N) {
NV thr = pts[i].score;
size_t j = i;
while (j < N && pts[j].score == thr) { if (pts[j].lab) TP += 1.0; else FP += 1.0; j++; }
NV sens = TP / P, spec = (Nn - FP) / Nn;
HV *restrict pt = newHV();
hv_stores(pt, "threshold", newSVnv(thr));
hv_stores(pt, "sensitivity", newSVnv(sens));
hv_stores(pt, "specificity", newSVnv(spec));
av_push(curve, newRV_noinc((SV *)pt));
NV jstat = sens + spec - 1.0;
if (jstat > best_j) { best_j = jstat; best_thr = thr; best_sens = sens; best_spec = spec; }
i = j;
}
Safefree(pts);
HV *restrict youden = newHV();
hv_stores(youden, "threshold", newSVnv(best_thr));
hv_stores(youden, "sensitivity", newSVnv(best_sens));
hv_stores(youden, "specificity", newSVnv(best_spec));
hv_stores(youden, "j", newSVnv(best_j));
HV *restrict ret = newHV();
hv_stores(ret, "auc", newSVnv(auc_val));
hv_stores(ret, "auc_se", newSVnv(se));
{ AV *restrict ci = newAV(); av_push(ci, newSVnv(lo)); av_push(ci, newSVnv(hi));
hv_stores(ret, "auc_ci", newRV_noinc((SV *)ci)); }
hv_stores(ret, "conf_level", newSVnv(conf_level));
hv_stores(ret, "n_pos", newSViv((IV)m));
hv_stores(ret, "n_neg", newSViv((IV)n));
hv_stores(ret, "n", newSViv((IV)N));
hv_stores(ret, "direction", newSVpv(lower_pos ? "<" : ">", 1));
hv_stores(ret, "youden", newRV_noinc((SV *)youden));
hv_stores(ret, "curve", newRV_noinc((SV *)curve));
hv_stores(ret, "method", newSVpv("ROC curve with DeLong AUC", 0));
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void bedroc(...)
PPCODE:
{
/*Boltzmann-Enhanced Discrimination of ROC (Truchon & Bayly 2007, eq. 36).
Rewards early recognition: actives ranked near the top count far more
than actives buried deep in the list. alpha sets how sharply the weight
decays with rank; ties get the average (mid)rank.*/
if (items == 1 && !SvROK(ST(0))) { //bedroc('h'|'H'|'?') => help
const char *restrict h = SvPV_nolen(ST(0));
if (strEQ(h, "h") || strEQ(h, "H") || strEQ(h, "?")) {
GV *ogv = gv_fetchpvs("STDOUT", 0, SVt_PVIO);
PerlIO *pio = (ogv && GvIO(ogv) && IoOFP(GvIO(ogv)))
? IoOFP(GvIO(ogv)) : PerlIO_stdout();
PerlIO_printf(pio,
"bedroc - Boltzmann-Enhanced Discrimination of ROC (Truchon & Bayly 2007)\n"
"\n"
"Early-recognition metric: scores actives ranked near the TOP of the list\n"
"far more than actives buried deep. Result is in [0, 1] (1 = ideal early\n"
"recognition, 0 = worst, ~0.5 = random-ish depending on alpha & R_a).\n"
"\n"
"USAGE\n"
" my $r = bedroc(\\@scores, \\@labels, alpha => 20);\n"
" print $r->{bedroc};\n"
"\n"
"ARGUMENTS\n"
" \\@scores ranking scores (higher = better by default)\n"
" \\@labels class labels, OR a numeric column when cutoff/active_frac\n"
" is given (then bedroc binarizes it for you)\n"
" alpha => 20 early-recognition weight, > 0 (Truchon-Bayly default 20)\n"
" positive => 1 label value that marks an active (string compare)\n"
" cutoff => x define actives as \\@labels entries with value >= x\n"
" active_frac=>f binarize \\@labels: take a fraction f in (0,1) as active,\n"
" from one tail (see active_side). Exactly ceil(f*N) actives,\n"
" clamped so both classes exist -- never dies for lack of a\n"
" pre-built 0/1 label. Mutually exclusive with cutoff.\n"
" active_side=> with active_frac: 'high' (default) = largest values are\n"
" 'high' active (like cutoff); 'low' = smallest values (e.g. actives\n"
" = strongest binders when the column is dG).\n"
" direction=>'>' '>' higher score ranks first (default); '<' flips\n"
" top => 0.05 also report enrichment in the top fraction (0..1]\n"
"\n"
"Unlike the common Python implementations (which need a pre-built 0/1 label\n"
"array, or reimplement a regression variant per script), active_frac lets one\n"
"call turn a raw measured column straight into a BEDROC score.\n"
"\n"
"RETURNS a hashref: bedroc, alpha, rie, rie_min, rie_max, n, n_active,\n"
" n_inactive, ra, direction, method, and (with top=>) enrichment =>\n"
" { fraction, n_top, active_count, expected, enrichment_factor }.\n"
"\n"
" bedroc('h'), bedroc('H') or bedroc('?') prints this help.\n");
XSRETURN(0);
}
}
if (items < 2 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV
|| !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV)
croak("Usage: bedroc(\\@scores, \\@labels, alpha => 20, "
"positive => 1, cutoff => x, active_frac => 0.1, "
"active_side => 'high', direction => '>', top => 0.05)");
NV alpha = 20.0; const char *restrict positive = "1"; int lower_pos = 0;
bool have_cutoff = 0, have_top = 0, have_frac = 0;
NV cutoff = 0.0, top = 0.0, active_frac = 0.0;
int frac_low = 0;
have_cutoff ? " (check cutoff)" : "");
}
qsort(pts, N, sizeof(ROCPt), rocpt_cmp_desc);
// sum over actives of exp(-alpha * midrank / N), 1-based ranks, best = 1 */
NV sum = 0.0;
for (size_t i = 0; i < N; ) {
size_t j = i;
while (j < N && pts[j].score == pts[i].score) j++;
NV midrank = ((NV)(i + 1) + (NV)j) / 2.0; // avg of positions i+1..j
NV w = exp(-alpha * midrank / (NV)N);
for (size_t k = i; k < j; k++) if (pts[k].lab) sum += w;
i = j;
}
NV ra = (NV)m / (NV)N;
NV rand_ = ra * (1.0 - exp(-alpha)) / (exp(alpha / (NV)N) - 1.0);
NV rie = sum / rand_;
NV f1 = ra * sinh(alpha / 2.0)
/ (cosh(alpha / 2.0) - cosh(alpha / 2.0 - alpha * ra));
NV f2 = 1.0 / (1.0 - exp(alpha * (1.0 - ra)));
NV bedroc = rie * f1 + f2;
NV rie_max = (1.0 - exp(-alpha * ra)) / (ra * (1.0 - exp(-alpha)));
NV rie_min = (1.0 - exp( alpha * ra)) / (ra * (1.0 - exp( alpha)));
HV *restrict ret = newHV();
hv_stores(ret, "bedroc", newSVnv(bedroc));
hv_stores(ret, "alpha", newSVnv(alpha));
hv_stores(ret, "rie", newSVnv(rie));
hv_stores(ret, "rie_min", newSVnv(rie_min));
hv_stores(ret, "rie_max", newSVnv(rie_max));
hv_stores(ret, "n", newSViv((IV)N));
hv_stores(ret, "n_active", newSViv((IV)m));
hv_stores(ret, "n_inactive", newSViv((IV)n));
hv_stores(ret, "ra", newSVnv(ra));
hv_stores(ret, "direction", newSVpv(lower_pos ? "<" : ">", 1));
hv_stores(ret, "method", newSVpv("BEDROC (Truchon-Bayly early recognition)", 0));
if (have_top) {//enrichment in the top fraction: EF = (hits/n_top) / R_a
size_t n_top = (size_t)ceil(top * (NV)N);
if (n_top < 1) n_top = 1; if (n_top > N) n_top = N;
size_t hits = 0;
for (size_t i = 0; i < n_top; i++) if (pts[i].lab) hits++;
NV expected = ra * (NV)n_top;
HV *enr = newHV();
hv_stores(enr, "fraction", newSVnv(top));
hv_stores(enr, "n_top", newSViv((IV)n_top));
hv_stores(enr, "active_count", newSViv((IV)hits));
hv_stores(enr, "expected", newSVnv(expected));
hv_stores(enr, "enrichment_factor",
newSVnv((expected > 0.0) ? ((NV)hits / (NV)n_top) / ra : NAN));
hv_stores(ret, "enrichment", newRV_noinc((SV *)enr));
}
Safefree(pts);
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void survfit(...)
PPCODE:
{
if (items < 2 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV
|| !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV)
croak("Usage: survfit(\\@time, \\@status, group => \\@grp, conf_level => 0.95)");
AV *restrict gav = NULL; NV conf_level = 0.95;
for (int i = 2; i + 1 < items; i += 2) {
const char *restrict k = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(k, "group")) {
if (!SvROK(v) || SvTYPE(SvRV(v)) != SVt_PVAV) croak("survfit: group must be an array ref");
gav = (AV *)SvRV(v);
}
else if (strEQ(k, "conf_level") || strEQ(k, "conf.level")) conf_level = SvNV(v);
else croak("survfit: unknown argument '%s'", k);
}
if (!(conf_level > 0.0 && conf_level < 1.0)) croak("survfit: conf_level must be between 0 and 1");
AV *labels = (AV *)sv_2mortal((SV *)newAV());
size_t N; SurvObs *o = srv_read(aTHX_ (AV *)SvRV(ST(0)), (AV *)SvRV(ST(1)), gav, &N, labels, "survfit");
qsort(o, N, sizeof(SurvObs), survobs_cmp);
SSize_t G = av_len(labels) + 1;
NV z = inverse_normal_cdf(1.0 - (1.0 - conf_level) / 2.0);
HV *strata = newHV();
for (SSize_t g = 0; g < G; g++) {// group size
size_t ng = 0; for (size_t i = 0; i < N; i++) if (o[i].grp == g) ng++;
AV *t_av=newAV(), *nr_av=newAV(), *ne_av=newAV(), *nc_av=newAV(),
*s_av=newAV(), *se_av=newAV(), *lo_av=newAV(), *hi_av=newAV();
NV S = 1.0, vterm = 0.0, median = NAN;
size_t at_risk = ng, total_events = 0;
size_t i = 0;
while (i < N) {
if (o[i].grp != g) { i++; continue; }
NV t = o[i].time;
size_t d = 0, c = 0, block = 0;
size_t j = i;
while (j < N && o[j].time == t) {
if (o[j].grp == g) { block++; if (o[j].status) d++; else c++; }
j++;
}
size_t nr = at_risk; // at risk just before t
if (d > 0 && nr > d) {
S *= 1.0 - (NV)d / (NV)nr;
vterm += (NV)d / ((NV)nr * (NV)(nr - d));
total_events += d;
} else if (d > 0) { // everyone remaining has an event
S = 0.0; total_events += d;
}
NV se_S = S * sqrt(vterm);
NV lo = (S > 0.0) ? S * exp(-z * sqrt(vterm)) : 0.0;
NV hi = (S > 0.0) ? S * exp( z * sqrt(vterm)) : 0.0;
if (hi > 1.0) hi = 1.0;
av_push(t_av, newSVnv(t));
av_push(nr_av, newSViv((IV)nr));
av_push(ne_av, newSViv((IV)d));
av_push(nc_av, newSViv((IV)c));
av_push(s_av, newSVnv(S));
av_push(se_av, newSVnv(se_S));
av_push(lo_av, newSVnv(lo));
av_push(hi_av, newSVnv(hi));
if (isnan(median) && S <= 0.5) median = t;
at_risk -= block;
i = j;
}
HV *st = newHV();
hv_stores(st, "time", newRV_noinc((SV *)t_av));
hv_stores(st, "n_risk", newRV_noinc((SV *)nr_av));
hv_stores(st, "n_event", newRV_noinc((SV *)ne_av));
hv_stores(st, "n_censor", newRV_noinc((SV *)nc_av));
hv_stores(st, "surv", newRV_noinc((SV *)s_av));
hv_stores(st, "std_err", newRV_noinc((SV *)se_av));
hv_stores(st, "lower", newRV_noinc((SV *)lo_av));
hv_stores(st, "upper", newRV_noinc((SV *)hi_av));
hv_stores(st, "n", newSViv((IV)ng));
hv_stores(st, "events", newSViv((IV)total_events));
hv_stores(st, "median", isnan(median) ? newSV(0) : newSVnv(median));
SV *key = *av_fetch(labels, g, 0);
STRLEN klen; const char *kp = SvPV(key, klen);
hv_store(strata, kp, klen, newRV_noinc((SV *)st), 0);
}
Safefree(o);
HV *ret = newHV();
hv_stores(ret, "strata", newRV_noinc((SV *)strata));
hv_stores(ret, "groups", newRV_inc((SV *)labels));
hv_stores(ret, "conf_level", newSVnv(conf_level));
hv_stores(ret, "method", newSVpv("Kaplan-Meier survival estimate", 0));
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void logrank_test(...)
PPCODE:
{
if (items < 3 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV
|| !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV
|| !SvROK(ST(2)) || SvTYPE(SvRV(ST(2))) != SVt_PVAV)
croak("Usage: logrank_test(\\@time, \\@status, \\@group)");
AV *restrict labels = (AV *)sv_2mortal((SV *)newAV());
size_t N; SurvObs *o = srv_read(aTHX_ (AV *)SvRV(ST(0)), (AV *)SvRV(ST(1)),
(AV *)SvRV(ST(2)), &N, labels, "logrank_test");
SSize_t G = av_len(labels) + 1;
if (G < 2) { Safefree(o); croak("logrank_test: need at least two groups"); }
qsort(o, N, sizeof(SurvObs), survobs_cmp);
NV *restrict O, *restrict E, *restrict V; Newx(O, G, NV); Newx(E, G, NV); Newx(V, G * G, NV);
for (SSize_t k = 0; k < G; k++) { O[k] = 0.0; E[k] = 0.0; }
for (SSize_t k = 0; k < G * G; k++) V[k] = 0.0;
size_t *restrict nrisk; Newx(nrisk, G, size_t);
for (SSize_t k = 0; k < G; k++) { size_t c = 0; for (size_t i = 0; i < N; i++) if (o[i].grp == k) c++; nrisk[k] = c; }
size_t i = 0;
while (i < N) {
NV t = o[i].time;
size_t d_tot = 0, n_tot = 0;
size_t *restrict dj; Newx(dj, G, size_t); for (SSize_t k = 0; k < G; k++) dj[k] = 0;
for (SSize_t k = 0; k < G; k++) n_tot += nrisk[k];
size_t j = i, block_j0 = 0; (void)block_j0;
while (j < N && o[j].time == t) { if (o[j].status) { dj[o[j].grp]++; d_tot++; } j++; }
if (d_tot > 0 && n_tot > 1) {
for (SSize_t a = 0; a < G; a++) {
NV nja = (NV)nrisk[a];
O[a] += (NV)dj[a];
E[a] += (NV)d_tot * nja / (NV)n_tot;
for (SSize_t b = 0; b < G; b++) {
NV njb = (NV)nrisk[b];
NV term = (NV)d_tot * ((NV)n_tot - (NV)d_tot) / ((NV)n_tot - 1.0)
* (((a == b) ? nja / (NV)n_tot : 0.0) - nja * njb / ((NV)n_tot * (NV)n_tot));
V[a * G + b] += term;
}
}
}
// remove this time's observations from the risk sets
for (size_t k2 = i; k2 < j; k2++) nrisk[o[k2].grp]--;
Safefree(dj);
i = j;
}
int m = (int)G - 1; // reduced dimension
NV *restrict Vr; Newx(Vr, m * m, NV);
NV *restrict OE; Newx(OE, m, NV);
for (int a = 0; a < m; a++) { OE[a] = O[a] - E[a]; for (int b = 0; b < m; b++) Vr[a*m+b] = V[a*G+b]; }
NV *restrict xsol; Newx(xsol, m, NV);
NV chi = 0.0;
if (srv_solve(Vr, OE, m, xsol) == 0)
for (int a = 0; a < m; a++) chi += OE[a] * xsol[a];
NV pval = get_p_value(chi, m);
HV *restrict ret = newHV();
hv_stores(ret, "statistic", newSVnv(chi));
hv_stores(ret, "parameter", newSViv(m));
hv_stores(ret, "p_value", newSVnv(pval));
AV *obs = newAV(), *exp_av = newAV();
for (SSize_t k = 0; k < G; k++) { av_push(obs, newSVnv(O[k])); av_push(exp_av, newSVnv(E[k])); }
hv_stores(ret, "observed", newRV_noinc((SV *)obs));
hv_stores(ret, "expected", newRV_noinc((SV *)exp_av));
hv_stores(ret, "groups", newRV_inc((SV *)labels));
hv_stores(ret, "method", newSVpv("Log-rank (Mantel-Cox) test", 0));
Safefree(o); Safefree(O); Safefree(E); Safefree(V); Safefree(nrisk);
Safefree(Vr); Safefree(OE); Safefree(xsol);
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void coxph(...)
PPCODE:
{
if (items < 3 || !SvROK(ST(0)) || SvTYPE(SvRV(ST(0))) != SVt_PVAV
|| !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV
|| !SvROK(ST(2)) || SvTYPE(SvRV(ST(2))) != SVt_PVAV)
croak("Usage: coxph(\\@time, \\@status, \\@covariate | [\\@x1, \\@x2, ...], "
"conf_level => 0.95, ties => 'efron', names => [...])");
AV *restrict tav = (AV *)SvRV(ST(0)), *sav = (AV *)SvRV(ST(1)), *Xav = (AV *)SvRV(ST(2));
SSize_t n = av_len(tav) + 1;
if (n < 2) croak("coxph: need at least two observations");
if (av_len(sav) + 1 != n) croak("coxph: time and status must be the same length");
//covariates: [\@x1, \@x2, ...] (multiple) or a single flat \@x
int p, multi = 0;
SV **restrict first = av_fetch(Xav, 0, 0);
if (first && *first && SvROK(*first) && SvTYPE(SvRV(*first)) == SVt_PVAV) { multi = 1; p = (int)(av_len(Xav) + 1); }
else { multi = 0; p = 1; }
if (p < 1) croak("coxph: need at least one covariate");
NV conf_level = 0.95; int breslow = 0, maxit = 25; NV eps = 1e-9;
AV *names_av = NULL;
for (int i = 3; i + 1 < items; i += 2) {
const char *k = SvPV_nolen(ST(i)); SV *v = ST(i + 1);
if (strEQ(k, "conf_level") || strEQ(k, "conf.level")) conf_level = SvNV(v);
else if (strEQ(k, "ties")) breslow = strEQ(SvPV_nolen(v), "breslow");
else if (strEQ(k, "maxit")) maxit = (int)SvIV(v);
else if (strEQ(k, "names")) { if (SvROK(v) && SvTYPE(SvRV(v)) == SVt_PVAV) names_av = (AV *)SvRV(v); }
else croak("coxph: unknown argument '%s'", k);
}
if (!(conf_level > 0.0 && conf_level < 1.0)) croak("coxph: conf_level must be between 0 and 1");
//pull data into contiguous C arrays
NV *restrict X; Newx(X, (size_t)n * p, NV);
NV *restrict tm; Newx(tm, n, NV);
int *restrict st; Newx(st, n, int);
for (SSize_t i = 0; i < n; i++) {
tm[i] = SvNV(*av_fetch(tav, i, 0));
st[i] = (SvNV(*av_fetch(sav, i, 0)) != 0.0) ? 1 : 0;
}
if (multi) {
for (int k = 0; k < p; k++) {
AV *col = (AV *)SvRV(*av_fetch(Xav, k, 0));
if (av_len(col) + 1 != n) { Safefree(X); Safefree(tm); Safefree(st); croak("coxph: covariate %d length mismatch", k + 1); }
for (SSize_t i = 0; i < n; i++) X[i * p + k] = SvNV(*av_fetch(col, i, 0));
}
} else {
if (av_len(Xav) + 1 != n) { Safefree(X); Safefree(tm); Safefree(st); croak("coxph: covariate length mismatch"); }
for (SSize_t i = 0; i < n; i++) X[i] = SvNV(*av_fetch(Xav, i, 0));
}
//observations sorted by time ascending (risk sets built time-descending)
TimeIdx *restrict ord; Newx(ord, n, TimeIdx);
for (SSize_t i = 0; i < n; i++) { ord[i].time = tm[i]; ord[i].idx = (int)i; }
qsort(ord, n, sizeof(TimeIdx), cmp_nv3);
NV *beta = NULL, *U = NULL, *Imat = NULL, *Iinv = NULL, *ratio = NULL,
*S1 = NULL, *S2 = NULL, *SD1 = NULL, *SD2 = NULL, *sumX_D = NULL, *eta = NULL, *w = NULL;
Newx(beta, p, NV); Newx(U, p, NV); Newx(Imat, p * p, NV); Newx(Iinv, p * p, NV);
Newx(ratio, p, NV); Newx(S1, p, NV); Newx(S2, p * p, NV);
Newx(SD1, p, NV); Newx(SD2, p * p, NV); Newx(sumX_D, p, NV);
Newx(eta, n, NV); Newx(w, n, NV);
//Newton update: beta += Iinv * U
for (int k = 0; k < p; k++) { NV s = 0.0; for (int l = 0; l < p; l++) s += Iinv[k * p + l] * U[l]; beta[k] += s; }
}
int nevent = 0; for (SSize_t i = 0; i < n; i++) if (st[i]) nevent++;
NV zc = inverse_normal_cdf(1.0 - (1.0 - conf_level) / 2.0);
AV *coef=newAV(), *hr=newAV(), *se=newAV(), *zv=newAV(), *pv=newAV(),
*ci=newAV(), *nm=newAV();
for (int k = 0; k < p; k++) {
NV b = beta[k];
NV sek = (Iinv[k * p + k] > 0.0) ? sqrt(Iinv[k * p + k]) : NAN;
NV zk = b / sek;
NV pk = 2.0 * approx_pnorm(-fabs(zk));
av_push(coef, newSVnv(b));
av_push(hr, newSVnv(exp(b)));
av_push(se, newSVnv(sek));
av_push(zv, newSVnv(zk));
av_push(pv, newSVnv(pk));
AV *cik = newAV();
av_push(cik, newSVnv(exp(b - zc * sek)));
av_push(cik, newSVnv(exp(b + zc * sek)));
av_push(ci, newRV_noinc((SV *)cik));
if (names_av && k <= av_len(names_av)) av_push(nm, newSVsv(*av_fetch(names_av, k, 0)));
else { SV *dn = newSVpvf("x%d", k + 1); av_push(nm, dn); }
}
NV lr = 2.0 * (loglik - loglik_null);
NV lr_p = get_p_value(lr, p);
HV *restrict ret = newHV();
hv_stores(ret, "coef", newRV_noinc((SV *)coef));
hv_stores(ret, "exp_coef", newRV_noinc((SV *)hr)); //hazard ratios
hv_stores(ret, "se", newRV_noinc((SV *)se));
hv_stores(ret, "z", newRV_noinc((SV *)zv));
hv_stores(ret, "p_value", newRV_noinc((SV *)pv));
hv_stores(ret, "conf_int", newRV_noinc((SV *)ci)); //on the HR scale
hv_stores(ret, "names", newRV_noinc((SV *)nm));
hv_stores(ret, "loglik", newSVnv(loglik));
hv_stores(ret, "loglik_null", newSVnv(loglik_null));
hv_stores(ret, "lr_stat", newSVnv(lr));
hv_stores(ret, "lr_df", newSViv(p));
hv_stores(ret, "lr_p_value", newSVnv(lr_p));
hv_stores(ret, "n", newSViv((IV)n));
hv_stores(ret, "nevent", newSViv(nevent));
hv_stores(ret, "iterations", newSViv(iter));
hv_stores(ret, "converged", newSViv(converged));
hv_stores(ret, "conf_level", newSVnv(conf_level));
hv_stores(ret, "ties", newSVpv(breslow ? "breslow" : "efron", 0));
hv_stores(ret, "method", newSVpv("Cox proportional hazards model", 0));
Safefree(X); Safefree(tm); Safefree(st); Safefree(ord);
Safefree(beta); Safefree(U); Safefree(Imat); Safefree(Iinv); Safefree(ratio);
Safefree(S1); Safefree(S2); Safefree(SD1); Safefree(SD2); Safefree(sumX_D);
Safefree(eta); Safefree(w);
ST(0) = sv_2mortal(newRV_noinc((SV *)ret));
XSRETURN(1);
}
void p_adjust(...)
PROTOTYPE: $;@
PPCODE:
if (items < 1)
croak("Usage: p_adjust($p_values, $method, columns => ...)");
SV *restrict p_sv = ST(0);
const char *restrict method = "holm";
SV *restrict cols_sv = NULL;
IV first_pair = 1;
/*The method may still arrive positionally, the way it always has;
anything after it (or after the frame) is key => value.*/
if (items > 1 && ((items - 1) % 2) == 1) {
if (!SvOK(ST(1)) || SvROK(ST(1)))
croak("p_adjust: the second argument must be an adjustment method name");
method = SvPV_nolen(ST(1));
first_pair = 2;
}
for (IV a = first_pair; a + 1 < (IV)items; a += 2) {
const char *restrict key = SvPV_nolen(ST(a));
SV *restrict val = ST(a + 1);
if (strEQ(key, "method")) method = SvPV_nolen(val);
else if (strEQ(key, "columns") || strEQ(key, "column")
|| strEQ(key, "cols") || strEQ(key, "col")) cols_sv = val;
else croak("p_adjust: unknown argument '%s'", key);
}
char meth[PA_METH_LEN];
pa_method(method, meth);
if (!pa_known(meth)) croak("Unknown p-value adjustment method: %s", method);
/*Which columns hold p-values? Nothing named means all of them. The
value is a flag, set once the column turns up in the frame.*/
HV *restrict want = NULL;
if (cols_sv && SvOK(cols_sv)) {
want = (HV*)sv_2mortal((SV*)newHV());
if (SvROK(cols_sv) && SvTYPE(SvRV(cols_sv)) == SVt_PVAV) {
AV *restrict cav = (AV*)SvRV(cols_sv);
for (SSize_t i = 0; i <= av_len(cav); i++) {
SV **restrict c = av_fetch(cav, i, 0);
if (!c || !SvOK(*c))
croak("p_adjust: undefined column name in 'columns'");
(void)hv_store_ent(want, *c, newSViv(0), 0);
}
} else if (SvROK(cols_sv)) {
croak("p_adjust: 'columns' must be a column name or an ARRAY "
"reference of column names");
} else {
(void)hv_store_ent(want, cols_sv, newSViv(0), 0);
}
if (HvUSEDKEYS(want) == 0) croak("p_adjust: 'columns' names no columns");
}
//Which of the five shapes is this?
enum { PA_FLAT, PA_AOA, PA_AOH, PA_HOA, PA_HOH } kind = PA_FLAT;
SV *restrict ref = SvROK(p_sv) ? SvRV(p_sv) : NULL;
if (!ref || (SvTYPE(ref) != SVt_PVAV && SvTYPE(ref) != SVt_PVHV))
croak("p_adjust: first argument must be an ARRAY reference of p-values, "
"or a reference to an AoA, AoH, HoA or HoH data frame");
if (SvTYPE(ref) == SVt_PVAV) {
AV *restrict av = (AV*)ref;
for (SSize_t i = 0; i <= av_len(av); i++) {
SV **restrict e = av_fetch(av, i, 0);
if (!e || !SvOK(*e)) continue; //undef p-value: still flat
below would read off a null pointer. av_fetch hands back
a deferred PVLV rather than the value, and SvOK on that is
false until its get-magic runs: without SvGETMAGIC every
element of a tied array looks undefined and this croaks
on data that is perfectly well defined.*/
for (size_t j = 0; j < len; j++) {
SV** restrict tv = av_fetch(av, j, 0);
if (tv) SvGETMAGIC(*tv);
if (tv && SvOK(*tv)) {
nums[k++] = SvNV(*tv);
} else {
if (nums != stackbuf) Safefree(nums);
/*UVuf, not %zu: croak() runs perl's own formatter, which does not
understand the C99 z modifier and prints it literally on older
perls (5.10 and 5.12 both do)*/
croak("median: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
} else {
/*AvARRAY, not av_fetch: the length is known and the cells
are right there, so the bounds check and the call per
element buy nothing*/
SV** restrict src = AvARRAY(av);
for (size_t j = 0; j < len; j++) {
SV* restrict tv = src[j];
if (tv && SvOK(tv)) {
nums[k++] = SvNV(tv);
} else {
if (nums != stackbuf) Safefree(nums);
croak("median: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
}
}
}
} else if (SvOK(arg)) {
nums[k++] = SvNV(arg);
} else {
if (nums != stackbuf) Safefree(nums);
croak("median: undefined value at argument index %" UVuf, (UV)i);
}
}
/*Select the middle value(s) rather than sorting all of them. For an
even count the lower of the pair is the largest value left below the
upper one, which a scan of that side finds without a second select.*/
if (total_count & 1) {
nv_select(nums, total_count, total_count / 2);
median_val = nums[total_count / 2];
} else {
const size_t up = total_count / 2;
nv_select(nums, total_count, up);
NV lower = nums[0];
for (size_t i = 1; i < up; i++) if (nums[i] > lower) lower = nums[i];
median_val = (lower + nums[up]) / 2.0;
}
if (nums != stackbuf) Safefree(nums);
RETVAL = median_val;
OUTPUT:
RETVAL
void intersection(...)
PROTOTYPE: @
PPCODE:
if (items == 0)
croak("intersection needs >= 1 array ref");
SP = set_multiplicity(aTHX_ SP, &ST(0), (size_t)items, 1, 0,
"intersection", GIMME_V);
SV* cor(SV* x_sv, SV* y_sv = &PL_sv_undef, const char* method = "pearson")
INIT:
// --- validate method
if (strcmp(method, "pearson") != 0 &&
strcmp(method, "spearman") != 0 &&
strcmp(method, "kendall") != 0)
croak("cor: unknown method '%s' (use 'pearson', 'spearman', or 'kendall')",
method);
// --- validate x
if (!SvROK(x_sv) || SvTYPE(SvRV(x_sv)) != SVt_PVAV)
croak("cor: x must be an ARRAY reference");
AV*restrict x_av = (AV*)SvRV(x_sv);
size_t nx = av_len(x_av) + 1;
if (nx == 0) croak("cor: x is empty");
// --- detect whether x is a flat vector or a matrix (AoA)
bool x_is_matrix = 0;
{
SV**restrict fp = av_fetch(x_av, 0, 0);
if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
x_is_matrix = 1;
}
// --- detect y
bool has_y = (SvOK(y_sv) && SvROK(y_sv) &&
SvTYPE(SvRV(y_sv)) == SVt_PVAV);
AV*restrict y_av = has_y ? (AV*)SvRV(y_sv) : NULL;
size_t ny = has_y ? av_len(y_av) + 1 : 0;
bool y_is_matrix = 0;
if (has_y && ny > 0) {
SV**restrict fp = av_fetch(y_av, 0, 0);
if (fp && SvROK(*fp) && SvTYPE(SvRV(*fp)) == SVt_PVAV)
y_is_matrix = 1;
}
CODE:
// Branch 1: both inputs are flat vectors → scalar result
if (!x_is_matrix && !y_is_matrix) {
if (!has_y) {
// cor(vector) == 1 by definition
RETVAL = newSVnv(1.0);
} else {
if (nx != ny)
croak("cor: x and y must have the same length (%lu vs %lu)",
nx, ny);
if (nx < 2)
croak("cor: need at least 2 observations");
NV *restrict xd, *restrict yd;
Newx(xd, nx, NV);
Newx(yd, ny, NV);
bool x_sd0 = 1, y_sd0 = 1;
symmetric = 1;
}
if (nrows < 2)
croak("cor: need at least 2 observations (got %lu)", nrows);
// -- build cache for symmetric case: compute upper triangle, store results, mirror to lower triangle
AV*restrict result_av = newAV();
av_extend(result_av, ncols_x - 1);
// Allocate per-row AVs up front so we can fill them in order
AV **restrict rows_out;
Newx(rows_out, ncols_x, AV*);
for (size_t i = 0; i < ncols_x; i++) {
rows_out[i] = newAV();
av_extend(rows_out[i], ncols_y - 1);
}
if (symmetric) {
// Upper triangle + diagonal, then mirror. r_cache[i][j] (j >= i) holds the computed value
NV **restrict r_cache;
Newx(r_cache, ncols_x, NV*);
for (size_t i = 0; i < ncols_x; i++)
Newx(r_cache[i], ncols_x, NV);
for (size_t i = 0; i < ncols_x; i++) {
r_cache[i][i] = 1.0; // diagonal
for (size_t j = i + 1; j < ncols_x; j++) {
NV r = compute_cor(col_x[i], col_x[j], nrows, method);
r_cache[i][j] = r;
r_cache[j][i] = r; // symmetry
}
}
// fill output AoA from cache
for (size_t i = 0; i < ncols_x; i++)
for (size_t j = 0; j < ncols_x; j++)
av_store(rows_out[i], j, newSVnv(r_cache[i][j]));
for (size_t i = 0; i < ncols_x; i++) Safefree(r_cache[i]);
Safefree(r_cache); r_cache = NULL;
} else {
// cross-correlation: every (i,j) pair is independent
for (size_t i = 0; i < ncols_x; i++)
for (size_t j = 0; j < ncols_y; j++)
av_store(rows_out[i], j, newSVnv(compute_cor(col_x[i], col_y[j], nrows, method)));
}
// push row AVs into result
for (size_t i = 0; i < ncols_x; i++)
av_store(result_av, i, newRV_noinc((SV*)rows_out[i]));
Safefree(rows_out); rows_out = NULL;
// -- free column arrays -------------------------------------
for (size_t j = 0; j < ncols_x; j++) Safefree(col_x[j]);
Safefree(col_x); col_x = NULL;
if (!symmetric) {
for (size_t j = 0; j < ncols_y; j++) Safefree(col_y[j]);
Safefree(col_y);
}
RETVAL = newRV_noinc((SV*)result_av);
}
OUTPUT:
RETVAL
void scale(...)
PROTOTYPE: @
PPCODE:
{
bool do_center_mean = TRUE, do_scale_sd = TRUE;
NV center_val = 0.0, scale_val = 1.0;
size_t data_items = items;
// 1. Parse Options Hash (if it exists as the last argument)
if (items > 0) {
SV*restrict last_arg = ST(items - 1);
if (SvROK(last_arg) && SvTYPE(SvRV(last_arg)) == SVt_PVHV) {
data_items = items - 1; // Exclude hash from data processing
HV*restrict opt_hv = (HV*)SvRV(last_arg);
// --- Parse 'center'
SV**restrict center_sv = hv_fetch(opt_hv, "center", 6, 0);
if (center_sv) {
SV*restrict val_sv = *center_sv;
if (!SvOK(val_sv)) {
do_center_mean = FALSE; center_val = 0.0;
} else {
char *restrict str = SvPV_nolen(val_sv);
//Trap booleans and empty strings before numeric checks
if (strcasecmp(str, "mean") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
do_center_mean = TRUE;
} else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
do_center_mean = FALSE; center_val = 0.0;
} else if (looks_like_number(val_sv)) {
do_center_mean = FALSE; center_val = SvNV(val_sv);
} else if (SvTRUE(val_sv)) {
do_center_mean = TRUE;
} else {
do_center_mean = FALSE; center_val = 0.0;
}
}
}
// --- Parse 'scale' ---
SV**restrict scale_sv = hv_fetch(opt_hv, "scale", 5, 0);
if (scale_sv) {
SV*restrict val_sv = *scale_sv;
if (!SvOK(val_sv)) {
do_scale_sd = FALSE; scale_val = 1.0;
} else {
char *restrict str = SvPV_nolen(val_sv);
if (strcasecmp(str, "sd") == 0 || strcasecmp(str, "true") == 0 || strcmp(str, "1") == 0) {
do_scale_sd = TRUE;
} else if (strcasecmp(str, "none") == 0 || strcasecmp(str, "false") == 0 || strcmp(str, "0") == 0 || strcmp(str, "") == 0) {
do_scale_sd = FALSE; scale_val = 1.0;
} else if (looks_like_number(val_sv)) {
do_scale_sd = FALSE; scale_val = SvNV(val_sv);
if (scale_val == 0.0) scale_val = 1.0; //Prevent Division By Zero
} else if (SvTRUE(val_sv)) {
do_scale_sd = TRUE;
} else {
do_scale_sd = FALSE; scale_val = 1.0;
}
}
}
}
}
// 2. Detect if the input is a Matrix (Array of Arrays)
bool is_matrix = FALSE;
if (data_items == 1) {
SV*restrict first_arg = ST(0);
}
} else if (final_rank == df_int) {
r_squared = 0.0; adj_r_squared = 0.0;
}
for (j = 0; j < p; j++) {
const char *restrict cname = design->col[j].name;
hv_store(coef_hv, cname, strlen(cname), newSVnv(beta[j]), 0);
av_push(terms_av, newSVpv(cname, 0));
HV *restrict row_hv = newHV();
if (aliased[j]) {
hv_store(row_hv, "Estimate", 8, newSVpv("NaN", 0), 0);
hv_store(row_hv, "Std. Error", 10, newSVpv("NaN", 0), 0);
hv_store(row_hv, "t value", 7, newSVpv("NaN", 0), 0);
hv_store(row_hv, "Pr(>|t|)", 8, newSVpv("NaN", 0), 0);
} else {
NV se = sqrt(rse_sq * XtX[j * p + j]);
NV t_val = (se > 0.0) ? (beta[j] / se) : (INFINITY * (beta[j] >= 0.0 ? 1.0 : -1.0));
NV p_val = get_t_pvalue(t_val, df_res, "two.sided");
hv_store(row_hv, "Estimate", 8, newSVnv(beta[j]), 0);
hv_store(row_hv, "Std. Error", 10, newSVnv(se), 0);
hv_store(row_hv, "t value", 7, newSVnv(t_val), 0);
hv_store(row_hv, "Pr(>|t|)", 8, newSVnv(p_val), 0);
}
hv_store(summary_hv, cname, strlen(cname), newRV_noinc((SV*)row_hv), 0);
}
hv_store(res_hv, "coefficients", 12, newRV_noinc((SV*)coef_hv), 0);
hv_store(res_hv, "fitted.values", 13, newRV_noinc((SV*)fitted_hv), 0);
hv_store(res_hv, "residuals", 9, newRV_noinc((SV*)resid_hv), 0);
hv_store(res_hv, "df.residual", 11, newSVuv(df_res), 0);
hv_store(res_hv, "rank", 4, newSVuv(final_rank), 0);
hv_store(res_hv, "rss", 3, newSVnv(rss), 0);
hv_store(res_hv, "summary", 7, newRV_noinc((SV*)summary_hv),0);
hv_store(res_hv, "terms", 5, newRV_noinc((SV*)terms_av), 0);
hv_store(res_hv, "r.squared", 9, newSVnv(r_squared), 0);
hv_store(res_hv, "adj.r.squared", 13, newSVnv(adj_r_squared), 0);
hv_store(res_hv, "xlevels", 7, newRV_inc((SV*)xlevels_hv), 0);
if (!isnan(f_stat)) {
AV *fstat_av = newAV();
av_push(fstat_av, newSVnv(f_stat));
av_push(fstat_av, newSViv(numdf));
av_push(fstat_av, newSViv(df_res));
hv_store(res_hv, "fstatistic", 10, newRV_noinc((SV*)fstat_av), 0);
hv_store(res_hv, "f.pvalue", 8, newSVnv(f_pvalue), 0);
}
for (i = 0; i < num_terms; i++) Safefree(terms[i]); Safefree(terms);
for (i = 0; i < num_uniq; i++) Safefree(uniq_terms[i]); Safefree(uniq_terms);
lm_design_free(aTHX_ design);
Safefree(X); Safefree(Y); Safefree(XtX); Safefree(XtY);
Safefree(beta); Safefree(aliased);
if (row_hashes) Safefree(row_hashes);
RETVAL = newRV_noinc((SV*)res_hv);
}
OUTPUT:
RETVAL
void seq(from, to, by = 1.0)
NV from
NV to
NV by
PPCODE:
{
if (by == 0.0) {//Handle the zero 'by' case
if (from == to) {
EXTEND(SP, 1);
mPUSHn(from);
XSRETURN(1);
} else {
croak("invalid 'by' argument: cannot be zero when from != to");
}
}
// Check for wrong direction / infinite loop
if ((from < to && by < 0.0) || (from > to && by > 0.0)) {
croak("wrong sign in 'by' argument");
}
/*Calculate number of elements.
R uses a small epsilon (like 1e-10) to avoid dropping the last
element due to floating point inaccuracies.*/
NV n_elements_d = (to - from) / by;
if (n_elements_d < 0.0) n_elements_d = 0.0;
size_t n_elements = (n_elements_d + 1e-10) + 1;
// Pre-extend the stack to avoid reallocating inside the loop
EXTEND(SP, n_elements);
for (size_t i = 0; i < n_elements; i++) {
mPUSHn(from + i * by);
}
XSRETURN(n_elements);
}
SV* rnorm(...)
CODE:
{
// Auto-seed the PRNG if the Perl script hasn't done so yet
AUTO_SEED_PRNG();
size_t n = 0;
NV mean = 0.0, sd = 1.0;
int arg_start = 0;
// Check if the first argument is a simple integer (rnorm(33))
if (items > 0 && SvIOK(ST(0)) && (items == 1 || items % 2 != 0)) {
n = (unsigned int)SvUV(ST(0));
arg_start = 1; // Start parsing named arguments from the second element
}
// --- Parse remaining named arguments from the flat stack ---
if ((items - arg_start) % 2 != 0) {
croak("Usage: rnorm(n), rnorm(n => 10, mean => 0, sd => 1), or rnorm(33, mean => 0)");
}
for (int i = arg_start; i < items; i += 2) {
const char* restrict key = SvPV_nolen(ST(i));
SV* restrict val = ST(i + 1);
if (strEQ(key, "n")) n = (unsigned int)SvUV(val);
else if (strEQ(key, "mean")) mean = SvNV(val);
else if (strEQ(key, "sd")) sd = SvNV(val);
else croak("rnorm: unknown argument '%s'", key);
}
if (sd < 0.0) croak("rnorm: standard deviation must be non-negative");
AV *restrict result_av = newAV();
if (n > 0) {
av_extend(result_av, n - 1);
Safefree(idx);
} else {
for (size_t i = 0; i < (size_t)n; i++)
av_push(ret_av, newSV(0));
}
ret = newRV_noinc((SV *)ret_av);
}
}
RETVAL = ret;
OUTPUT:
RETVAL
SV* dnorm(...)
CODE:
{
if (items < 1) {
croak("Usage: dnorm(x), dnorm(x, mean => 0, sd => 1, log => 0)");
}
SV*restrict x_sv = ST(0);
NV mean = 0.0, sd = 1.0; //defaults
bool give_log = 0;
// --- Parse remaining named arguments from the flat stack ---
if ((items - 1) % 2 != 0) {
croak("dnorm: Expected an even number of key-value named arguments after 'x'");
}
for (size_t i = 1; i < items; i += 2) {
const char* restrict key = SvPV_nolen(ST(i));
SV* restrict val = ST(i + 1);
if (strEQ(key, "mean")) mean = SvNV(val);
else if (strEQ(key, "sd")) sd = SvNV(val);
else if (strEQ(key, "log")) give_log = SvTRUE(val) ? 1 : 0;
else croak("dnorm: unknown argument '%s'", key);
}
// --- Branch based on scalar vs. arrayref for 'x' ---
if (SvROK(x_sv) && SvTYPE(SvRV(x_sv)) == SVt_PVAV) {
// x is an array reference
AV *restrict x_av = (AV*)SvRV(x_sv);
IV n = av_len(x_av) + 1;
AV *restrict result_av = newAV();
if (n > 0) {
av_extend(result_av, n - 1);
for (IV i = 0; i < n; i++) {
SV **restrict elem = av_fetch(x_av, i, 0);
NV x_val = (elem && *elem) ? SvNV(*elem) : NAN;
NV res = c_dnorm(x_val, mean, sd, give_log);
av_store(result_av, i, newSVnv(res));
}
}
RETVAL = newRV_noinc((SV*)result_av);
} else {
// x is a single numeric scalar
NV x_val = SvNV(x_sv);
NV res = c_dnorm(x_val, mean, sd, give_log);
RETVAL = newSVnv(res);
}
}
OUTPUT:
RETVAL
void merge(...)
PPCODE:
{
if (items < 2)
croak("Usage: merge($left, $right, how => 'inner'|'left'|'right'|"
"'outer'|'cross', on => 'col' | ['c1','c2'] "
"[, 'left.on' => .., 'right.on' => ..] "
"[, suffixes => ['.x','.y']] [, 'output.type' => 'aoh'|'hoa'])");
if ((items - 2) & 1)
croak("merge: options after the two frames must be name => value pairs");
SV *restrict left = ST(0);
SV *restrict right = ST(1);
SV *restrict how_sv = NULL, *restrict on_sv = NULL;
SV *restrict lon_sv = NULL, *restrict ron_sv = NULL;
SV *restrict suf_sv = NULL, *restrict out_sv = NULL;
for (int oi = 2; oi < items; oi += 2) {
STRLEN ol;
const char *restrict on = SvPV(ST(oi), ol);
SV *restrict ov = ST(oi + 1);
if (strEQ(on, "how")) how_sv = ov;
else if (strEQ(on, "on") || strEQ(on, "by")) on_sv = ov;
else if (strEQ(on, "left.on") || strEQ(on, "left_on")
|| strEQ(on, "by.x")) lon_sv = ov;
else if (strEQ(on, "right.on") || strEQ(on, "right_on")
|| strEQ(on, "by.y")) ron_sv = ov;
else if (strEQ(on, "suffixes")) suf_sv = ov;
else if (strEQ(on, "output.type") || strEQ(on, "output_type")
|| strEQ(on, "out")) out_sv = ov;
else croak("merge: unknown option '%s'", on);
}
//how
int how = MG_INNER;
if (how_sv && SvOK(how_sv)) {
const char *restrict h = SvPV_nolen(how_sv);
if (strEQ(h, "inner")) how = MG_INNER;
else if (strEQ(h, "left")) how = MG_LEFT;
else if (strEQ(h, "right")) how = MG_RIGHT;
else if (strEQ(h, "outer") || strEQ(h, "full")) how = MG_OUTER;
else if (strEQ(h, "cross")) how = MG_CROSS;
else croak("merge: how must be 'inner', 'left', 'right', 'outer', or "
"'cross' (got '%s')", h);
}
if (on_sv && (lon_sv || ron_sv))
croak("merge: give either 'on'/'by' or 'left.on'/'right.on', not both");
if ((lon_sv && !ron_sv) || (ron_sv && !lon_sv))
croak("merge: 'left.on' and 'right.on' must be given together");
if (how == MG_CROSS && (on_sv || lon_sv || ron_sv))
croak("merge: a cross join takes no join keys");
ENTER; SAVETMPS;
//suffixes
SV *restrict suf0 = NULL, *restrict suf1 = NULL;
if (suf_sv) {
if (!SvROK(suf_sv) || SvTYPE(SvRV(suf_sv)) != SVt_PVAV
|| av_len((AV *)SvRV(suf_sv)) != 1)
croak("merge: suffixes must be a two-element array-ref, e.g. ['.x','.y']");
AV *restrict sa = (AV *)SvRV(suf_sv);
Newx(kv, ncols ? ncols : 1, SV *);
SAVEFREEPV(kv);
Newx(cv, ncols ? ncols : 1, AV *);
SAVEFREEPV(cv);
//one pass to collect columns and find the longest
n = 0;
ci = 0;
hv_iterinit(in);
while ((he = hv_iternext(in))) {
SV *restrict val = HeVAL(he);
size_t len;
if (!val || !SvROK(val) || SvTYPE(SvRV(val)) != SVt_PVAV)
croak("hoa2hoh: column '%s' is not an arrayref",
SvPV_nolen(hv_iterkeysv(he)));
kv[ci] = hv_iterkeysv(he); //mortal; valid until our LEAVE
cv[ci] = (AV *)SvRV(val);
len = (size_t)(av_len(cv[ci]) + 1);
if (len > n)
n = len;
ci++;
}
ncols = ci;
out = newHV();
sv_2mortal((SV *)out); //reclaimed on croak; +1'd below on success
for (i = 0; i < n; i++) {
HV *restrict row;
SV *restrict rowname;
SV **restrict kp = av_fetch(keycol, i, 0);
if (!kp || !*kp || !SvOK(*kp))
croak("hoa2hoh: key column '%s' has an undefined value at row %" UVuf,
SvPV_nolen(key), (UV)i);
rowname = *kp;
if (hv_exists_ent(out, rowname, 0))
croak("hoa2hoh: duplicate row name '%s'", SvPV_nolen(rowname));
row = newHV();
for (ci = 0; ci < ncols; ci++) {
SV **restrict cp = av_fetch(cv[ci], i, 0);
SV *restrict cell = (cp && *cp) ? newSVsv(*cp) : newSV(0);
(void)hv_store_ent(row, kv[ci], cell, 0);
}
(void)hv_store_ent(out, rowname, newRV_noinc((SV *)row), 0);
}
RETVAL = newRV_inc((SV *)out);
FREETMPS;
LEAVE;
}
OUTPUT:
RETVAL
void vals(data, colname_sv)
SV *data
SV *colname_sv
PREINIT:
bool is_aoh = 0, is_hoh = 0;
const char *restrict colname = NULL;
STRLEN collen = 0;
AV *restrict src_av = NULL;
HV *restrict src_hv = NULL;
SSize_t n = 0;
AV *restrict out_av = NULL;
PPCODE:
{
if (!SvOK(colname_sv))
croak("vals: column name must be defined");
colname = SvPV(colname_sv, collen); //kept for the error message
if (!SvROK(data))
croak("vals: first argument must be an array-ref (AoH) or hash-ref (HoA, HoH)");
//---- classify $data: AoH (arrayref) vs HoA/HoH (hashref) --------
if (SvTYPE(SvRV(data)) == SVt_PVAV) {
is_aoh = 1;
src_av = (AV *)SvRV(data);
n = av_len(src_av) + 1;
} else if (SvTYPE(SvRV(data)) == SVt_PVHV) {
src_hv = (HV *)SvRV(data);
hv_iterinit(src_hv);
HE *restrict he = hv_iternext(src_hv);
if (he) {
SV *restrict val = HeVAL(he);
if (val && SvROK(val) && SvTYPE(SvRV(val)) == SVt_PVHV)
is_hoh = 1; //a hash whose values are hashes => HoH
//else leave is_aoh/is_hoh = 0 => HoA path below
}
// empty hash: is_aoh = is_hoh = 0 => HoA path yields []
} else {
croak("vals: first argument must be an array-ref (AoH) or hash-ref (HoA, HoH)");
}
//out_av is mortalised up front so any later croak frees it cleanly
out_av = newAV();
sv_2mortal((SV *)out_av);
if (is_aoh) { // AoH
if (n > 0) av_extend(out_av, n - 1);
for (SSize_t i = 0; i < n; i++) {
SV **restrict rp = av_fetch(src_av, i, 0);
SV *restrict row = (rp && *rp) ? *rp : &PL_sv_undef;
/* strict: a row must be a hash-ref, else fail here with the index
rather than returning undef and letting the caller die vaguely*/
if (!SvOK(row))
croak("vals: AoH row %" IVdf " is undef (expected a hash-ref)", (IV)i);
if (!SvROK(row) || SvTYPE(SvRV(row)) != SVt_PVHV)
croak("vals: AoH row %" IVdf " is not a hash-ref", (IV)i);
HE *restrict ent = hv_fetch_ent((HV *)SvRV(row), colname_sv, 0, 0);
// a valid row that simply lacks the column still yields undef (R-like NA)
SV *restrict cell = (ent && HeVAL(ent)) ? HeVAL(ent) : &PL_sv_undef;
/*copy, so the result is independent of the source and undef
slots are writable (not the shared read-only PL_sv_undef)*/
av_push(out_av, newSVsv(cell));
}
} else if (is_hoh) { // HoH
n = hv_iterinit(src_hv);
if (n > 0) {
av_extend(out_av, n - 1);
ENTER;
SV **restrict keys; SV **restrict rows;
Newx(keys, n, SV *); SAVEFREEPV(keys);
Newx(rows, n, SV *); SAVEFREEPV(rows);
SSize_t cnt = 0;
HE *restrict he;
while ((he = hv_iternext(src_hv)) && cnt < n) {
keys[cnt] = hv_iterkeysv(he); //mortal copy of the key
rows[cnt] = HeVAL(he);
keys[j + 1] = k;
rows[j + 1] = r;
}
for (SSize_t i = 0; i < cnt; i++) {
SV *restrict row_sv = rows[i];
// strict: name the offending key instead of silently emitting undef
if (!row_sv || !SvROK(row_sv) || SvTYPE(SvRV(row_sv)) != SVt_PVHV)
croak("vals: HoH value for key '%s' is not a hash-ref",
SvPV_nolen(keys[i]));
HE *restrict ent = hv_fetch_ent((HV *)SvRV(row_sv), colname_sv, 0, 0);
SV *restrict cell = (ent && HeVAL(ent)) ? HeVAL(ent) : &PL_sv_undef;
av_push(out_av, newSVsv(cell));
}
LEAVE;
}
} else { // HoA
if (hv_iterinit(src_hv) > 0) { //non-empty hash
HE *restrict colent = hv_fetch_ent(src_hv, colname_sv, 0, 0);
SV *restrict cv = colent ? HeVAL(colent) : NULL;
if (!cv || !SvROK(cv) || SvTYPE(SvRV(cv)) != SVt_PVAV)
croak("vals: column '%s' not found or is not an array-ref", colname);
AV *restrict col_av = (AV *)SvRV(cv);
n = av_len(col_av) + 1;
if (n > 0) {
/*the length is known, so the result is sized once and filled
straight through AvARRAY rather than pushed a cell at a time*/
av_extend(out_av, n - 1);
SV **restrict d = AvARRAY(out_av);
SV **restrict src = AvARRAY(col_av);
const SSize_t srcn = AvFILLp(col_av) + 1;
for (SSize_t i = 0; i < n; i++)
d[i] = newSVsv((i < srcn && src[i]) ? src[i] : &PL_sv_undef);
AvFILLp(out_av) = n - 1;
}
}
}
/*out_av is mortal (freed on any croak); newRV_inc balances that so the
returned RV holds the surviving reference -- newRV_noinc here would
double-free with the mortal.*/
XPUSHs(sv_2mortal(newRV_inc((SV *)out_av)));
XSRETURN(1);
}
void
_qcut_core(data_ref, probs_ref, drop_dups, want_codes)
SV *data_ref
SV *probs_ref
IV drop_dups
IV want_codes
PREINIT:
AV *data_av;
AV *probs_av;
AV *edge_av;
AV *code_av = NULL;
SV **el;
IV n, m, i, j, ne, w;
NV *srt = NULL;
NV *edges = NULL;
NV p, h, frac, v;
IV lo, bin, lo2, hi2, mid, k;
PPCODE:
if (!SvROK(data_ref) || SvTYPE(SvRV(data_ref)) != SVt_PVAV)
croak("_qcut_core: data must be an ARRAY reference");
if (!SvROK(probs_ref) || SvTYPE(SvRV(probs_ref)) != SVt_PVAV)
croak("_qcut_core: probs must be an ARRAY reference");
data_av = (AV *) SvRV(data_ref);
probs_av = (AV *) SvRV(probs_ref);
n = av_len(data_av) + 1;
m = av_len(probs_av) + 1;
if (n < 1)
croak("_qcut_core: need at least one data value");
if (m < 2)
croak("_qcut_core: need at least two probabilities (one bin)");
Newx(srt, n, NV);
for (i = 0; i < n; i++) {
el = av_fetch(data_av, i, 0);
srt[i] = (el && SvOK(*el)) ? SvNV(*el) : 0.0;
}
qsort(srt, (size_t) n, sizeof(NV), cmp_nv3);
//quantile cutpoints via linear interpolation (numpy/pandas default)
Newx(edges, m, NV);
for (j = 0; j < m; j++) {
el = av_fetch(probs_av, j, 0);
p = el ? SvNV(*el) : 0.0;
if (p < 0.0) p = 0.0;
if (p > 1.0) p = 1.0;
h = (NV)(n - 1) * p;
lo = (IV) floor((double) h);
frac = h - (NV) lo;
if (lo + 1 < n)
edges[j] = srt[lo] + frac * (srt[lo + 1] - srt[lo]);
else
edges[j] = srt[lo];
}
//guard fp drift: enforce non-decreasing edges
for (j = 1; j < m; j++)
if (edges[j] < edges[j - 1])
edges[j] = edges[j - 1];
Safefree(srt); //no longer needed once cutpoints exist
//duplicate edges: raise (default) or drop
w = 1;
for (j = 1; j < m; j++) {
if (edges[j] == edges[w - 1]) {
if (!drop_dups) {
Safefree(edges);
croak("_qcut_core: bin edges are not unique; pass duplicates => 'drop' (or use fewer bins)");
}
} else {
edges[w++] = edges[j];
}
}
ne = w;
if (ne < 2) {
Safefree(edges);
croak("_qcut_core: data has too few distinct values to form bins");
}
edge_av = newAV();
av_extend(edge_av, ne - 1);
for (j = 0; j < ne; j++)
av_push(edge_av, newSVnv(edges[j]));
/*assign each original value to a 0-based bin only if codes are wanted;
lowest bin is inclusive on both ends*/
if (want_codes) {
code_av = newAV();
av_extend(code_av, n - 1);
for (i = 0; i < n; i++) {
el = av_fetch(data_av, i, 0);
v = (el && SvOK(*el)) ? SvNV(*el) : 0.0;
if (v <= edges[0]) {
bin = 0;
} else if (v >= edges[ne - 1]) {
bin = ne - 2;
} else {
lo2 = 1;
hi2 = ne - 1;
k = ne - 1;
while (lo2 <= hi2) {
mid = lo2 + ((hi2 - lo2) >> 1);
if (edges[mid] >= v) {
k = mid;
hi2 = mid - 1;
} else {
lo2 = mid + 1;
}
}
bin = k - 1;
}
av_push(code_av, newSViv(bin));
}
}
Safefree(edges);
EXTEND(SP, 2);
if (want_codes)
PUSHs(sv_2mortal(newRV_noinc((SV *) code_av)));
else
PUSHs(&PL_sv_undef);
PUSHs(sv_2mortal(newRV_noinc((SV *) edge_av)));
void get_union(...)
PROTOTYPE: @
PREINIT:
HV*restrict seen;
AV*restrict order;
size_t nrefs, n, oi, olen;
int gimme;
PPCODE:
gimme = GIMME_V;
nrefs = items;
if (nrefs == 0)
croak("union needs >= 1 array ref");
seen = (HV*)sv_2mortal((SV*)newHV());
order = (AV*)sv_2mortal((SV*)newAV()); //buffer: pushing to the stack while still reading ST() would clobber the args
n = 0;
for (size_t i = 0; i < nrefs; i++) {
SV*restrict arg = ST(i);
AV*restrict av;
size_t len;
if (!(SvROK(arg) && SvTYPE(SvRV(arg)) == SVt_PVAV))
croak("union: argument index %" UVuf " of %" UVuf " total (max index %" UVuf ") is not an array reference", (UV)i, (UV)nrefs, (UV)(nrefs - 1));
av = (AV*)SvRV(arg);
len = (size_t)(av_len(av) + 1);
for (size_t j = 0; j < len; j++) {
SV**restrict tv = av_fetch(av, j, 0);
STRLEN klen;
const char*restrict key;
I32 hklen;
if (!(tv && SvOK(*tv)))
croak("union: undefined value at array ref index %" UVuf " (argument %" UVuf ")", (UV)j, (UV)i);
key = SvPV(*tv, klen);
hklen = SvUTF8(*tv) ? -(I32)klen : (I32)klen;
if (hv_exists(seen, key, hklen))
continue;
(void)hv_store(seen, key, hklen, &PL_sv_undef, 0);
n++;
if (gimme != G_SCALAR)
av_push(order, newSVsv(*tv));
}
}
if (gimme == G_SCALAR) {
XPUSHs(sv_2mortal(newSVuv(n)));
} else {
olen = (size_t)(av_len(order) + 1);
for (oi = 0; oi < olen; oi++) {
SV**restrict e = av_fetch(order, oi, 0);
if (e && *e)
XPUSHs(sv_2mortal(newSVsv(*e)));
}
}
void Lonly(...)
PROTOTYPE: @
PPCODE:
if (items == 0)
croak("Lonly needs >= 1 array ref");
SP = set_multiplicity(aTHX_ SP, &ST(0), (size_t)items, 0, 0,
"Lonly", GIMME_V);
void Ronly(...)
PROTOTYPE: @
PPCODE:
if (items == 0)
croak("Ronly needs >= 1 array ref");
/*mirror of Lonly: values only in the LAST array (from_last = 1), so
the two-array Ronly(a,b) still equals Lonly(b,a).*/
SP = set_multiplicity(aTHX_ SP, &ST(0), (size_t)items, 0, 1,
"Ronly", GIMME_V);
void is_equivalent(...)
PROTOTYPE: @
PPCODE:
if (items < 2)
croak("is_equivalent needs >= 2 array refs (got %" UVuf ")", (UV)items);
XPUSHs(sv_2mortal(newSViv(set_equivalent(aTHX_ &ST(0), (size_t)items, "is_equivalent"))));
SV* pnorm(...)
CODE:
{
if (items < 1)
croak("Usage: pnorm(x), pnorm(x, mean => 0, sd => 1, lower => 1, log => 0)");
SV *restrict x_sv = ST(0);
NV mean = 0.0, sd = 1.0; // defaults
bool lower_tail = 1, give_log = 0;
if ((items - 1) % 2 != 0)
croak("pnorm: Expected an even number of key-value named arguments after 'x'");
for (size_t i = 1; i < items; i += 2) {
const char *restrict key = SvPV_nolen(ST(i));
SV *restrict val = ST(i + 1);
if (strEQ(key, "mean")) mean = SvNV(val);
else if (strEQ(key, "sd")) sd = SvNV(val);
else if (strEQ(key, "lower") || strEQ(key, "lower.tail")) lower_tail = SvTRUE(val) ? 1 : 0;
else if (strEQ(key, "log") || strEQ(key, "log.p")) give_log = SvTRUE(val) ? 1 : 0;
else croak("pnorm: unknown argument '%s'", key);
}
if (sd < 0.0)
warn("pnorm: standard deviation must be non-negative");
if (SvROK(x_sv) && SvTYPE(SvRV(x_sv)) == SVt_PVAV) {
AV *restrict x_av = (AV*)SvRV(x_sv);
IV n = av_len(x_av) + 1;
AV *restrict result_av = newAV();
if (n > 0) {
av_extend(result_av, n - 1);
for (IV i = 0; i < n; i++) {
SV **restrict elem = av_fetch(x_av, i, 0);
NV x_val = (elem && *elem) ? SvNV(*elem) : NAN;
NV res = (NV)c_pnorm((double)x_val, (double)mean, (double)sd,
lower_tail, give_log);
av_store(result_av, i, newSVnv(res));
}
}
RETVAL = newRV_noinc((SV*)result_av);
} else {
NV x_val = SvNV(x_sv);
NV res = (NV)c_pnorm((double)x_val, (double)mean, (double)sd,
lower_tail, give_log);
RETVAL = newSVnv(res);
}
}
OUTPUT:
RETVAL
# Private numeric helpers. These replace pure-Perl ports that used to live in
# LikeR.pm (_lgamma/_igamc/_pchisq_upper); igamc() here is the one authoritative
# implementation, so the Perl and C copies can no longer drift apart. Not
# exported -- callers inside Stats::LikeR use them unqualified.
NV _igamc(a, x)
NV a
NV x
CODE:
RETVAL = igamc(a, x);
( run in 1.500 second using v1.01-cache-2.11-cpan-4e7a2411597 )