Mirrors/miller
1
0
Fork 0
mirror of https://github.com/johnkerl/miller.git synced 2026-07-19 01:15:21 +00:00
Code Issues Projects Releases Wiki Activity
miller/c/mapping/mapper_stats2.c

1009 lines
40 KiB
C
Raw Permalink Blame History

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "lib/mlrutil.h"
#include "lib/mlr_globals.h"
#include "lib/mlrmath.h"
#include "lib/mlrstat.h"
#include "containers/sllv.h"
#include "containers/slls.h"
#include "lib/string_array.h"
#include "containers/lhmslv.h"
#include "containers/lhms2v.h"
#include "containers/lhmsv.h"
#include "containers/mixutil.h"
#include "containers/dvector.h"
#include "mapping/mappers.h"
#include "cli/argparse.h"
typedef enum _bivar_measure_t {
DO_CORR,
DO_COV,
DO_COVX,
DO_LINREG_PCA
} bivar_measure_t;
// ----------------------------------------------------------------
struct _stats2_acc_t; // forward reference for method definitions
typedef void stats2_ingest_func_t(void* pvstate, double x, double y);
typedef void stats2_emit_func_t(void* pvstate, char* name1, char* name2, lrec_t* poutrec);
typedef void stats2_fit_func_t(void* pvstate, double x, double y, lrec_t* poutrec);
typedef void stats2_free_func_t(struct _stats2_acc_t* pstats2_acc);
typedef struct _stats2_acc_t {
void* pvstate;
stats2_ingest_func_t* pingest_func;
stats2_emit_func_t* pemit_func;
stats2_fit_func_t* pfit_func;
stats2_free_func_t* pfree_func; // virtual destructor
} stats2_acc_t;
typedef struct _mapper_stats2_state_t {
ap_state_t* pargp;
slls_t* paccumulator_names;
string_array_t* pvalue_field_name_pairs;
slls_t* pgroup_by_field_names;
lhmslv_t* acc_groups;
lhmslv_t* record_groups;
int do_verbose;
int do_iterative_stats;
int do_hold_and_fit;
} mapper_stats2_state_t;
typedef stats2_acc_t* stats2_alloc_func_t(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
// ----------------------------------------------------------------
static void mapper_stats2_usage(FILE* o, char* argv0, char* verb);
static mapper_t* mapper_stats2_parse_cli(int* pargi, int argc, char** argv,
cli_reader_opts_t* _, cli_writer_opts_t* __);
static mapper_t* mapper_stats2_alloc(ap_state_t* pargp, slls_t* paccumulator_names,
string_array_t* pvalue_field_name_pairs, slls_t* pgroup_by_field_names,
int do_verbose, int do_iterative_stats, int do_hold_and_fit);
static void mapper_stats2_free(mapper_t* pmapper, context_t* _);
static sllv_t* mapper_stats2_process(lrec_t* pinrec, context_t* pctx, void* pvstate);
static void mapper_stats2_ingest(lrec_t* pinrec, context_t* pctx, mapper_stats2_state_t* pstate);
static sllv_t* mapper_stats2_emit_all(mapper_stats2_state_t* pstate);
static void mapper_stats2_emit(mapper_stats2_state_t* pstate, lrec_t* pinrec,
char* value_field_name_1, char* value_field_name_2, lhmsv_t* pacc_fields_to_acc_state);
static sllv_t* mapper_stats2_fit_all(mapper_stats2_state_t* pstate);
static stats2_acc_t* make_stats2 (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_linreg_pca_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_linreg_ols_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_r2_alloc (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_logireg_alloc (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_corr_cov_alloc (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, bivar_measure_t do_which, int do_verbose);
static stats2_acc_t* stats2_corr_alloc (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_cov_alloc (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
static stats2_acc_t* stats2_covx_alloc (char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose);
// ----------------------------------------------------------------
typedef struct _stats2_acc_lookup_t {
char* name;
stats2_alloc_func_t* palloc_func;
char* desc;
} stats2_acc_lookup_t;
static stats2_acc_lookup_t stats2_acc_lookup_table[] = {
{"linreg-pca", stats2_linreg_pca_alloc, "Linear regression using principal component analysis"},
{"linreg-ols", stats2_linreg_ols_alloc, "Linear regression using ordinary least squares"},
{"r2", stats2_r2_alloc, "Quality metric for linreg-ols (linreg-pca emits its own)"},
{"logireg", stats2_logireg_alloc, "Logistic regression"},
{"corr", stats2_corr_alloc, "Sample correlation"},
{"cov", stats2_cov_alloc, "Sample covariance"},
{"covx", stats2_covx_alloc, "Sample-covariance matrix"},
};
static int stats2_acc_lookup_table_length = sizeof(stats2_acc_lookup_table) / sizeof(stats2_acc_lookup_table[0]);
// ----------------------------------------------------------------
mapper_setup_t mapper_stats2_setup = {
.verb = "stats2",
.pusage_func = mapper_stats2_usage,
.pparse_func = mapper_stats2_parse_cli,
.ignores_input = FALSE,
};
// ----------------------------------------------------------------
static void mapper_stats2_usage(FILE* o, char* argv0, char* verb) {
fprintf(o, "Usage: %s %s [options]\n", argv0, verb);
fprintf(o, "Computes bivariate statistics for one or more given field-name pairs,\n");
fprintf(o, "accumulated across the input record stream.\n");
fprintf(o, "-a {linreg-ols,corr,...} Names of accumulators: one or more of:\n");
for (int i = 0; i < stats2_acc_lookup_table_length; i++) {
fprintf(o, " %-12s %s\n", stats2_acc_lookup_table[i].name, stats2_acc_lookup_table[i].desc);
}
fprintf(o, "-f {a,b,c,d} Value-field name-pairs on which to compute statistics.\n");
fprintf(o, " There must be an even number of names.\n");
fprintf(o, "-g {e,f,g} Optional group-by-field names.\n");
fprintf(o, "-v Print additional output for linreg-pca.\n");
fprintf(o, "-s Print iterative stats. Useful in tail -f contexts (in which\n");
fprintf(o, " case please avoid pprint-format output since end of input\n");
fprintf(o, " stream will never be seen).\n");
fprintf(o, "--fit Rather than printing regression parameters, applies them to\n");
fprintf(o, " the input data to compute new fit fields. All input records are\n");
fprintf(o, " held in memory until end of input stream. Has effect only for\n");
fprintf(o, " linreg-ols, linreg-pca, and logireg.\n");
fprintf(o, "Only one of -s or --fit may be used.\n");
fprintf(o, "Example: %s %s -a linreg-pca -f x,y\n", argv0, verb);
fprintf(o, "Example: %s %s -a linreg-ols,r2 -f x,y -g size,shape\n", argv0, verb);
fprintf(o, "Example: %s %s -a corr -f x,y\n", argv0, verb);
}
static mapper_t* mapper_stats2_parse_cli(int* pargi, int argc, char** argv,
cli_reader_opts_t* _, cli_writer_opts_t* __)
{
slls_t* paccumulator_names = NULL;
string_array_t* pvalue_field_names = NULL;
slls_t* pgroup_by_field_names = slls_alloc();
int do_verbose = FALSE;
int do_iterative_stats = FALSE;
int do_hold_and_fit = FALSE;
int allow_int_float = TRUE;
char* verb = argv[(*pargi)++];
ap_state_t* pstate = ap_alloc();
ap_define_string_list_flag(pstate, "-a", &paccumulator_names);
ap_define_string_array_flag(pstate, "-f", &pvalue_field_names);
ap_define_string_list_flag(pstate, "-g", &pgroup_by_field_names);
ap_define_true_flag(pstate, "-v", &do_verbose);
ap_define_true_flag(pstate, "-s", &do_iterative_stats);
ap_define_true_flag(pstate, "--fit", &do_hold_and_fit);
// The -F isn't used for stats2: all arithmetic here is floating-point. Yet
// it is supported for step and stats1 for all applicable stats1/step
// accumulators, so we accept here as well for all applicable stats2
// accumulators (i.e. none of them).
ap_define_false_flag(pstate, "-F", &allow_int_float);
if (!ap_parse(pstate, verb, pargi, argc, argv)) {
mapper_stats2_usage(stderr, argv[0], verb);
return NULL;
}
if (do_iterative_stats && do_hold_and_fit) {
mapper_stats2_usage(stderr, argv[0], verb);
return NULL;
}
if (paccumulator_names == NULL || pvalue_field_names == NULL) {
mapper_stats2_usage(stderr, argv[0], verb);
return NULL;
}
if ((pvalue_field_names->length % 2) != 0) {
mapper_stats2_usage(stderr, argv[0], verb);
return NULL;
}
return mapper_stats2_alloc(pstate, paccumulator_names, pvalue_field_names, pgroup_by_field_names,
do_verbose, do_iterative_stats, do_hold_and_fit);
}
// ----------------------------------------------------------------
static mapper_t* mapper_stats2_alloc(ap_state_t* pargp, slls_t* paccumulator_names,
string_array_t* pvalue_field_name_pairs, slls_t* pgroup_by_field_names,
int do_verbose, int do_iterative_stats, int do_hold_and_fit)
{
mapper_t* pmapper = mlr_malloc_or_die(sizeof(mapper_t));
mapper_stats2_state_t* pstate = mlr_malloc_or_die(sizeof(mapper_stats2_state_t));
pstate->pargp = pargp;
pstate->paccumulator_names = paccumulator_names;
pstate->pvalue_field_name_pairs = pvalue_field_name_pairs; // caller validates length is even
pstate->pgroup_by_field_names = pgroup_by_field_names;
pstate->acc_groups = lhmslv_alloc();
pstate->record_groups = lhmslv_alloc();
pstate->do_verbose = do_verbose;
pstate->do_iterative_stats = do_iterative_stats;
pstate->do_hold_and_fit = do_hold_and_fit;
pmapper->pvstate = pstate;
pmapper->pprocess_func = mapper_stats2_process;
pmapper->pfree_func = mapper_stats2_free;
return pmapper;
}
static void mapper_stats2_free(mapper_t* pmapper, context_t* _) {
mapper_stats2_state_t* pstate = pmapper->pvstate;
slls_free(pstate->paccumulator_names);
string_array_free(pstate->pvalue_field_name_pairs);
slls_free(pstate->pgroup_by_field_names);
// lhmslv_free and lhmsv_free will free the hashmap keys; we need to free
// the void-star hashmap values.
for (lhmslve_t* pa = pstate->acc_groups->phead; pa != NULL; pa = pa->pnext) {
lhms2v_t* pgroup_to_acc_field = pa->pvvalue;
for (lhms2ve_t* pb = pgroup_to_acc_field->phead; pb != NULL; pb = pb->pnext) {
lhmsv_t* pacc_fields_to_acc_state = pb->pvvalue;
for (lhmsve_t* pc = pacc_fields_to_acc_state->phead; pc != NULL; pc = pc->pnext) {
stats2_acc_t* pstats2_acc = pc->pvvalue;
pstats2_acc->pfree_func(pstats2_acc);
}
lhmsv_free(pacc_fields_to_acc_state);
}
lhms2v_free(pgroup_to_acc_field);
}
lhmslv_free(pstate->acc_groups);
for (lhmslve_t* pd = pstate->record_groups->phead; pd != NULL; pd = pd->pnext) {
sllv_t* plist = pd->pvvalue;
sllv_free(plist);
}
lhmslv_free(pstate->record_groups);
ap_free(pstate->pargp);
free(pstate);
free(pmapper);
}
// ================================================================
// Given: accumulate corr,cov on values x,y group by a,b.
// Example input: Example output:
// a b x y a b x_corr x_cov y_corr y_cov
// s t 1 2 s t 2 6 2 8
// u v 3 4 u v 1 3 1 4
// s t 5 6 u w 1 7 1 9
// u w 7 9
//
// Multilevel hashmap structure:
// {
// ["s","t"] : { <--- group-by field names
// ["x","y"] : { <--- value field names
// "corr" : stats2_corr_t object,
// "cov" : stats2_cov_t object
// }
// },
// ["u","v"] : {
// ["x","y"] : {
// "corr" : stats2_corr_t object,
// "cov" : stats2_cov_t object
// }
// },
// ["u","w"] : {
// ["x","y"] : {
// "corr" : stats2_corr_t object,
// "cov" : stats2_cov_t object
// }
// },
// }
// ================================================================
// In the iterative case, add to the current record its current group's stats fields.
// In the non-iterative case, produce output only at the end of the input stream.
static sllv_t* mapper_stats2_process(lrec_t* pinrec, context_t* pctx, void* pvstate) {
mapper_stats2_state_t* pstate = pvstate;
if (pinrec != NULL) {
mapper_stats2_ingest(pinrec, pctx, pstate);
if (pstate->do_iterative_stats) {
// The input record is modified in this case, with new fields appended
return sllv_single(pinrec);
} else if (pstate->do_hold_and_fit) {
// The input record is held by the ingestor
return NULL;
} else {
lrec_free(pinrec);
return NULL;
}
} else if (!pstate->do_iterative_stats) {
if (!pstate->do_hold_and_fit) {
return mapper_stats2_emit_all(pstate);
} else {
return mapper_stats2_fit_all(pstate);
}
} else {
return NULL;
}
}
// ----------------------------------------------------------------
static void mapper_stats2_ingest(lrec_t* pinrec, context_t* pctx, mapper_stats2_state_t* pstate) {
// ["s", "t"]
slls_t* pgroup_by_field_values = mlr_reference_selected_values_from_record(pinrec, pstate->pgroup_by_field_names);
if (pgroup_by_field_values == NULL) {
return;
}
lhms2v_t* pgroup_to_acc_field = lhmslv_get(pstate->acc_groups, pgroup_by_field_values);
if (pgroup_to_acc_field == NULL) {
pgroup_to_acc_field = lhms2v_alloc();
lhmslv_put(pstate->acc_groups, slls_copy(pgroup_by_field_values), pgroup_to_acc_field, FREE_ENTRY_KEY);
}
if (pstate->do_hold_and_fit) { // Retain the input record in memory, for fitting and delivery at end of stream
sllv_t* group_to_records = lhmslv_get(pstate->record_groups, pgroup_by_field_values);
if (group_to_records == NULL) {
group_to_records = sllv_alloc();
lhmslv_put(pstate->record_groups, slls_copy(pgroup_by_field_values), group_to_records, FREE_ENTRY_KEY);
}
sllv_append(group_to_records, pinrec);
}
// for [["x","y"]]
int n = pstate->pvalue_field_name_pairs->length;
for (int i = 0; i < n; i += 2) {
char* value_field_name_1 = pstate->pvalue_field_name_pairs->strings[i];
char* value_field_name_2 = pstate->pvalue_field_name_pairs->strings[i+1];
lhmsv_t* pacc_fields_to_acc_state = lhms2v_get(pgroup_to_acc_field, value_field_name_1, value_field_name_2);
if (pacc_fields_to_acc_state == NULL) {
pacc_fields_to_acc_state = lhmsv_alloc();
lhms2v_put(pgroup_to_acc_field, value_field_name_1, value_field_name_2, pacc_fields_to_acc_state, NO_FREE);
}
char* sval1 = lrec_get(pinrec, value_field_name_1);
char* sval2 = lrec_get(pinrec, value_field_name_2);
if (sval1 == NULL) // Key not present
continue;
if (*sval1 == 0) // Key present with null value
continue;
if (sval2 == NULL) // Key not present
continue;
if (*sval2 == 0) // Key present with null value
continue;
// for ["corr", "cov"]
sllse_t* pc = pstate->paccumulator_names->phead;
for ( ; pc != NULL; pc = pc->pnext) {
char* stats2_acc_name = pc->value;
stats2_acc_t* pstats2_acc = lhmsv_get(pacc_fields_to_acc_state, stats2_acc_name);
if (pstats2_acc == NULL) {
pstats2_acc = make_stats2(value_field_name_1, value_field_name_2, stats2_acc_name, pstate->do_verbose);
if (pstats2_acc == NULL) {
fprintf(stderr, "mlr stats2: accumulator \"%s\" not found.\n",
stats2_acc_name);
exit(1);
}
lhmsv_put(pacc_fields_to_acc_state, stats2_acc_name, pstats2_acc, NO_FREE);
}
if (sval1 == NULL || sval2 == NULL)
continue;
double dval1 = mlr_double_from_string_or_die(sval1);
double dval2 = mlr_double_from_string_or_die(sval2);
pstats2_acc->pingest_func(pstats2_acc->pvstate, dval1, dval2);
}
if (pstate->do_iterative_stats) {
mapper_stats2_emit(pstate, pinrec, value_field_name_1, value_field_name_2,
pacc_fields_to_acc_state);
}
}
slls_free(pgroup_by_field_values);
}
// ----------------------------------------------------------------
static sllv_t* mapper_stats2_emit_all(mapper_stats2_state_t* pstate) {
sllv_t* poutrecs = sllv_alloc();
for (lhmslve_t* pa = pstate->acc_groups->phead; pa != NULL; pa = pa->pnext) {
lrec_t* poutrec = lrec_unbacked_alloc();
// Add in a=s,b=t fields:
slls_t* pgroup_by_field_values = pa->key;
sllse_t* pb = pstate->pgroup_by_field_names->phead;
sllse_t* pc = pgroup_by_field_values->phead;
for ( ; pb != NULL && pc != NULL; pb = pb->pnext, pc = pc->pnext) {
lrec_put(poutrec, pb->value, pc->value, 0);
}
// Add in fields such as x_y_corr, etc.
lhms2v_t* pgroup_to_acc_field = pa->pvvalue;
// For "x","y"
for (lhms2ve_t* pd = pgroup_to_acc_field->phead; pd != NULL; pd = pd->pnext) {
char* value_field_name_1 = pd->key1;
char* value_field_name_2 = pd->key2;
lhmsv_t* pacc_fields_to_acc_state = pd->pvvalue;
mapper_stats2_emit(pstate, poutrec, value_field_name_1, value_field_name_2,
pacc_fields_to_acc_state);
// For "corr", "linreg"
for (lhmsve_t* pe = pacc_fields_to_acc_state->phead; pe != NULL; pe = pe->pnext) {
stats2_acc_t* pstats2_acc = pe->pvvalue;
pstats2_acc->pemit_func(pstats2_acc->pvstate, value_field_name_1, value_field_name_2, poutrec);
}
}
sllv_append(poutrecs, poutrec);
}
sllv_append(poutrecs, NULL);
return poutrecs;
}
static void mapper_stats2_emit(mapper_stats2_state_t* pstate, lrec_t* poutrec,
char* value_field_name_1, char* value_field_name_2, lhmsv_t* pacc_fields_to_acc_state)
{
// For "corr", "linreg"
for (lhmsve_t* pe = pacc_fields_to_acc_state->phead; pe != NULL; pe = pe->pnext) {
stats2_acc_t* pstats2_acc = pe->pvvalue;
pstats2_acc->pemit_func(pstats2_acc->pvstate, value_field_name_1, value_field_name_2, poutrec);
}
}
// ----------------------------------------------------------------
static sllv_t* mapper_stats2_fit_all(mapper_stats2_state_t* pstate) {
sllv_t* poutrecs = sllv_alloc();
for (lhmslve_t* pa = pstate->acc_groups->phead; pa != NULL; pa = pa->pnext) {
slls_t* pgroup_by_field_values = pa->key;
sllv_t* precords = lhmslv_get(pstate->record_groups, pgroup_by_field_values);
while (precords->phead) {
lrec_t* prec = sllv_pop(precords);
lhms2v_t* pgroup_to_acc_field = pa->pvvalue;
// For "x","y"
for (lhms2ve_t* pd = pgroup_to_acc_field->phead; pd != NULL; pd = pd->pnext) {
char* value_field_name_1 = pd->key1;
char* value_field_name_2 = pd->key2;
lhmsv_t* pacc_fields_to_acc_state = pd->pvvalue;
// For "linreg-ols", "logireg"
for (lhmsve_t* pe = pacc_fields_to_acc_state->phead; pe != NULL; pe = pe->pnext) {
stats2_acc_t* pstats2_acc = pe->pvvalue;
if (pstats2_acc->pfit_func != NULL) {
char* sx = lrec_get(prec, value_field_name_1);
char* sy = lrec_get(prec, value_field_name_2);
if (sx != NULL && sy != NULL) {
double x = mlr_double_from_string_or_die(sx);
double y = mlr_double_from_string_or_die(sy);
pstats2_acc->pfit_func(pstats2_acc->pvstate, x, y, prec);
}
}
}
}
sllv_append(poutrecs, prec);
}
}
sllv_append(poutrecs, NULL);
return poutrecs;
}
// ================================================================
// Given: accumulate corr,cov on values x,y group by a,b.
// Example input: Example output:
// a b x y a b x_corr x_cov y_corr y_cov
// s t 1 2 s t 2 6 2 8
// u v 3 4 u v 1 3 1 4
// s t 5 6 u w 1 7 1 9
// u w 7 9
//
// Multilevel hashmap structure:
// {
// ["s","t"] : { <--- group-by field names
// ["x","y"] : { <--- value field names
// "corr" : stats2_corr_t object,
// "cov" : stats2_cov_t object
// }
// },
// ["u","v"] : {
// ["x","y"] : {
// "corr" : stats2_corr_t object,
// "cov" : stats2_cov_t object
// }
// },
// ["u","w"] : {
// ["x","y"] : {
// "corr" : stats2_corr_t object,
// "cov" : stats2_cov_t object
// }
// },
// }
// ================================================================
// ----------------------------------------------------------------
static stats2_acc_t* make_stats2(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
for (int i = 0; i < stats2_acc_lookup_table_length; i++)
if (streq(stats2_acc_name, stats2_acc_lookup_table[i].name))
return stats2_acc_lookup_table[i].palloc_func(value_field_name_1, value_field_name_2, stats2_acc_name, do_verbose);
return NULL;
}
// ----------------------------------------------------------------
typedef struct _stats2_linreg_ols_state_t {
unsigned long long count;
double sumx;
double sumy;
double sumx2;
double sumxy;
char* m_output_field_name;
char* b_output_field_name;
char* n_output_field_name;
char* fit_output_field_name;
int fit_ready;
double m;
double b;
} stats2_linreg_ols_state_t;
static void stats2_linreg_ols_ingest(void* pvstate, double x, double y) {
stats2_linreg_ols_state_t* pstate = pvstate;
pstate->count++;
pstate->sumx += x;
pstate->sumy += y;
pstate->sumx2 += x*x;
pstate->sumxy += x*y;
}
static void stats2_linreg_ols_emit(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
stats2_linreg_ols_state_t* pstate = pvstate;
if (pstate->count < 2) {
lrec_put(poutrec, pstate->m_output_field_name, "", NO_FREE);
lrec_put(poutrec, pstate->b_output_field_name, "", NO_FREE);
} else {
double m, b;
mlr_get_linear_regression_ols(pstate->count, pstate->sumx, pstate->sumx2, pstate->sumxy, pstate->sumy, &m, &b);
char* mval = mlr_alloc_string_from_double(m, MLR_GLOBALS.ofmt);
char* bval = mlr_alloc_string_from_double(b, MLR_GLOBALS.ofmt);
lrec_put(poutrec, pstate->m_output_field_name, mval, FREE_ENTRY_VALUE);
lrec_put(poutrec, pstate->b_output_field_name, bval, FREE_ENTRY_VALUE);
}
char* nval = mlr_alloc_string_from_ll(pstate->count);
lrec_put(poutrec, pstate->n_output_field_name, nval, FREE_ENTRY_VALUE);
}
static void stats2_linreg_ols_fit(void* pvstate, double x, double y, lrec_t* poutrec) {
stats2_linreg_ols_state_t* pstate = pvstate;
if (!pstate->fit_ready) {
mlr_get_linear_regression_ols(pstate->count, pstate->sumx, pstate->sumx2, pstate->sumxy, pstate->sumy,
&pstate->m, &pstate->b);
pstate->fit_ready = TRUE;
}
if (pstate->count < 2) {
lrec_put(poutrec, pstate->fit_output_field_name, "", NO_FREE);
} else {
double yfit = pstate->m * x + pstate->b;
char* sfit = mlr_alloc_string_from_double(yfit, MLR_GLOBALS.ofmt);
lrec_put(poutrec, pstate->fit_output_field_name, sfit, FREE_ENTRY_VALUE);
}
}
static void stats2_linreg_ols_free(stats2_acc_t* pstats2_acc) {
stats2_linreg_ols_state_t* pstate = pstats2_acc->pvstate;
free(pstate->m_output_field_name);
free(pstate->b_output_field_name);
free(pstate->n_output_field_name);
free(pstate->fit_output_field_name);
free(pstate);
free(pstats2_acc);
}
static stats2_acc_t* stats2_linreg_ols_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
stats2_acc_t* pstats2_acc = mlr_malloc_or_die(sizeof(stats2_acc_t));
stats2_linreg_ols_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_linreg_ols_state_t));
pstate->count = 0LL;
pstate->sumx = 0.0;
pstate->sumy = 0.0;
pstate->sumx2 = 0.0;
pstate->sumxy = 0.0;
pstate->m_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_ols_m");
pstate->b_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_ols_b");
pstate->n_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_ols_n");
pstate->fit_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_ols_fit");
pstate->fit_ready = FALSE;
pstate->m = -999.0;
pstate->b = -999.0;
pstats2_acc->pvstate = (void*)pstate;
pstats2_acc->pingest_func = stats2_linreg_ols_ingest;
pstats2_acc->pemit_func = stats2_linreg_ols_emit;
pstats2_acc->pfit_func = stats2_linreg_ols_fit;
pstats2_acc->pfree_func = stats2_linreg_ols_free;
return pstats2_acc;
}
// ----------------------------------------------------------------
#define LOGIREG_DVECTOR_INITIAL_SIZE 1024
typedef struct _stats2_logireg_state_t {
dvector_t* pxs;
dvector_t* pys;
char* m_output_field_name;
char* b_output_field_name;
char* n_output_field_name;
char* fit_output_field_name;
int fit_ready;
double m;
double b;
} stats2_logireg_state_t;
static void stats2_logireg_ingest(void* pvstate, double x, double y) {
stats2_logireg_state_t* pstate = pvstate;
dvector_append(pstate->pxs, x);
dvector_append(pstate->pys, y);
}
static void stats2_logireg_emit(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
stats2_logireg_state_t* pstate = pvstate;
if (pstate->pxs->size < 2) {
lrec_put(poutrec, pstate->m_output_field_name, "", NO_FREE);
lrec_put(poutrec, pstate->b_output_field_name, "", NO_FREE);
} else {
double m, b;
mlr_logistic_regression(pstate->pxs->data, pstate->pys->data, pstate->pxs->size, &m, &b);
char* mval = mlr_alloc_string_from_double(m, MLR_GLOBALS.ofmt);
char* bval = mlr_alloc_string_from_double(b, MLR_GLOBALS.ofmt);
lrec_put(poutrec, pstate->m_output_field_name, mval, FREE_ENTRY_VALUE);
lrec_put(poutrec, pstate->b_output_field_name, bval, FREE_ENTRY_VALUE);
}
char* nval = mlr_alloc_string_from_ll(pstate->pxs->size);
lrec_put(poutrec, pstate->n_output_field_name, nval, FREE_ENTRY_VALUE);
}
static void stats2_logireg_free(stats2_acc_t* pstats2_acc) {
stats2_logireg_state_t* pstate = pstats2_acc->pvstate;
free(pstate->m_output_field_name);
free(pstate->b_output_field_name);
free(pstate->n_output_field_name);
free(pstate->fit_output_field_name);
dvector_free(pstate->pxs);
dvector_free(pstate->pys);
free(pstate);
free(pstats2_acc);
}
static void stats2_logireg_fit(void* pvstate, double x, double y, lrec_t* poutrec) {
stats2_logireg_state_t* pstate = pvstate;
if (!pstate->fit_ready) {
mlr_logistic_regression(pstate->pxs->data, pstate->pys->data, pstate->pxs->size, &pstate->m, &pstate->b);
pstate->fit_ready = TRUE;
}
if (pstate->pxs->size < 2) {
lrec_put(poutrec, pstate->fit_output_field_name, "", NO_FREE);
} else {
double yfit = 1.0 / (1.0 + exp(-pstate->m*x - pstate->b));
char* fitval = mlr_alloc_string_from_double(yfit, MLR_GLOBALS.ofmt);
lrec_put(poutrec, pstate->fit_output_field_name, fitval, FREE_ENTRY_VALUE);
}
}
static stats2_acc_t* stats2_logireg_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
stats2_acc_t* pstats2_acc = mlr_malloc_or_die(sizeof(stats2_acc_t));
stats2_logireg_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_logireg_state_t));
pstate->pxs = dvector_alloc(LOGIREG_DVECTOR_INITIAL_SIZE);
pstate->pys = dvector_alloc(LOGIREG_DVECTOR_INITIAL_SIZE);
pstate->m_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_logistic_m");
pstate->b_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_logistic_b");
pstate->n_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_logistic_n");
pstate->fit_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_logistic_fit");
pstate->fit_ready = FALSE;
pstate->m = -999.0;
pstate->b = -999.0;
pstats2_acc->pvstate = (void*)pstate;
pstats2_acc->pingest_func = stats2_logireg_ingest;
pstats2_acc->pemit_func = stats2_logireg_emit;
pstats2_acc->pfit_func = stats2_logireg_fit;
pstats2_acc->pfree_func = stats2_logireg_free;
return pstats2_acc;
}
// ----------------------------------------------------------------
// http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient
// Alternatively, just use sqrt(corr) as defined above.
typedef struct _stats2_r2_state_t {
unsigned long long count;
double sumx;
double sumy;
double sumx2;
double sumxy;
double sumy2;
char* r2_output_field_name;
} stats2_r2_state_t;
static void stats2_r2_ingest(void* pvstate, double x, double y) {
stats2_r2_state_t* pstate = pvstate;
pstate->count++;
pstate->sumx += x;
pstate->sumy += y;
pstate->sumx2 += x*x;
pstate->sumxy += x*y;
pstate->sumy2 += y*y;
}
static void stats2_r2_emit(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
stats2_r2_state_t* pstate = pvstate;
if (pstate->count < 2LL) {
lrec_put(poutrec, pstate->r2_output_field_name, "", NO_FREE);
} else {
unsigned long long n = pstate->count;
double sumx = pstate->sumx;
double sumy = pstate->sumy;
double sumx2 = pstate->sumx2;
double sumy2 = pstate->sumy2;
double sumxy = pstate->sumxy;
double numerator = n*sumxy - sumx*sumy;
numerator = numerator * numerator;
double denominator = (n*sumx2 - sumx*sumx) * (n*sumy2 - sumy*sumy);
double output = numerator/denominator;
char* val = mlr_alloc_string_from_double(output, MLR_GLOBALS.ofmt);
lrec_put(poutrec, pstate->r2_output_field_name, val, FREE_ENTRY_VALUE);
}
}
static void stats2_r2_free(stats2_acc_t* pstats2_acc) {
stats2_r2_state_t* pstate = pstats2_acc->pvstate;
free(pstate->r2_output_field_name);
free(pstate);
free(pstats2_acc);
}
static stats2_acc_t* stats2_r2_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
stats2_acc_t* pstats2_acc = mlr_malloc_or_die(sizeof(stats2_acc_t));
stats2_r2_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_r2_state_t));
pstate->count = 0LL;
pstate->sumx = 0.0;
pstate->sumy = 0.0;
pstate->sumx2 = 0.0;
pstate->sumxy = 0.0;
pstate->sumy2 = 0.0;
pstate->r2_output_field_name = mlr_paste_4_strings(value_field_name_1, "_", value_field_name_2, "_r2");
pstats2_acc->pvstate = (void*)pstate;
pstats2_acc->pingest_func = stats2_r2_ingest;
pstats2_acc->pemit_func = stats2_r2_emit;
pstats2_acc->pfit_func = NULL;
pstats2_acc->pfree_func = stats2_r2_free;
return pstats2_acc;
}
// ----------------------------------------------------------------
// Corr(X,Y) = Cov(X,Y) / sigma_X sigma_Y.
typedef struct _stats2_corr_cov_state_t {
unsigned long long count;
double sumx;
double sumy;
double sumx2;
double sumxy;
double sumy2;
bivar_measure_t do_which;
int do_verbose;
char* covx_00_output_field_name;
char* covx_01_output_field_name;
char* covx_10_output_field_name;
char* covx_11_output_field_name;
char* pca_m_output_field_name;
char* pca_b_output_field_name;
char* pca_n_output_field_name;
char* pca_q_output_field_name;
char* pca_l1_output_field_name;
char* pca_l2_output_field_name;
char* pca_v11_output_field_name;
char* pca_v12_output_field_name;
char* pca_v21_output_field_name;
char* pca_v22_output_field_name;
char* pca_fit_output_field_name;
int fit_ready;
double m;
double b;
double q;
char* corr_output_field_name;
char* cov_output_field_name;
} stats2_corr_cov_state_t;
static void stats2_corr_cov_ingest(void* pvstate, double x, double y) {
stats2_corr_cov_state_t* pstate = pvstate;
pstate->count++;
pstate->sumx += x;
pstate->sumy += y;
pstate->sumx2 += x*x;
pstate->sumxy += x*y;
pstate->sumy2 += y*y;
}
static void stats2_corr_cov_emit(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
stats2_corr_cov_state_t* pstate = pvstate;
if (pstate->do_which == DO_COVX) {
char* key00 = pstate->covx_00_output_field_name;
char* key01 = pstate->covx_01_output_field_name;
char* key10 = pstate->covx_10_output_field_name;
char* key11 = pstate->covx_11_output_field_name;
if (pstate->count < 2LL) {
lrec_put(poutrec, key00, "", NO_FREE);
lrec_put(poutrec, key01, "", NO_FREE);
lrec_put(poutrec, key10, "", NO_FREE);
lrec_put(poutrec, key11, "", NO_FREE);
} else {
double Q[2][2];
mlr_get_cov_matrix(pstate->count,
pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy, Q);
char* val00 = mlr_alloc_string_from_double(Q[0][0], MLR_GLOBALS.ofmt);
char* val01 = mlr_alloc_string_from_double(Q[0][1], MLR_GLOBALS.ofmt);
char* val10 = mlr_alloc_string_from_double(Q[1][0], MLR_GLOBALS.ofmt);
char* val11 = mlr_alloc_string_from_double(Q[1][1], MLR_GLOBALS.ofmt);
lrec_put(poutrec, key00, val00, FREE_ENTRY_VALUE);
lrec_put(poutrec, key01, val01, FREE_ENTRY_VALUE);
lrec_put(poutrec, key10, val10, FREE_ENTRY_VALUE);
lrec_put(poutrec, key11, val11, FREE_ENTRY_VALUE);
}
} else if (pstate->do_which == DO_LINREG_PCA) {
char* keym = pstate->pca_m_output_field_name;
char* keyb = pstate->pca_b_output_field_name;
char* keyn = pstate->pca_n_output_field_name;
char* keyq = pstate->pca_q_output_field_name;
char* keyl1 = pstate->pca_l1_output_field_name;
char* keyl2 = pstate->pca_l2_output_field_name;
char* keyv11 = pstate->pca_v11_output_field_name;
char* keyv12 = pstate->pca_v12_output_field_name;
char* keyv21 = pstate->pca_v21_output_field_name;
char* keyv22 = pstate->pca_v22_output_field_name;
if (pstate->count < 2LL) {
lrec_put(poutrec, keym, "", NO_FREE);
lrec_put(poutrec, keyb, "", NO_FREE);
lrec_put(poutrec, keyn, "", NO_FREE);
lrec_put(poutrec, keyq, "", NO_FREE);
if (pstate->do_verbose) {
lrec_put(poutrec, keyl1, "", NO_FREE);
lrec_put(poutrec, keyl2, "", NO_FREE);
lrec_put(poutrec, keyv11, "", NO_FREE);
lrec_put(poutrec, keyv12, "", NO_FREE);
lrec_put(poutrec, keyv21, "", NO_FREE);
lrec_put(poutrec, keyv22, "", NO_FREE);
}
} else {
double Q[2][2];
mlr_get_cov_matrix(pstate->count,
pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy, Q);
double l1, l2; // Eigenvalues
double v1[2], v2[2]; // Eigenvectors
mlr_get_real_symmetric_eigensystem(Q, &l1, &l2, v1, v2);
double x_mean = pstate->sumx / pstate->count;
double y_mean = pstate->sumy / pstate->count;
double m, b, q;
mlr_get_linear_regression_pca(l1, l2, v1, v2, x_mean, y_mean, &m, &b, &q);
lrec_put(poutrec, keym, mlr_alloc_string_from_double(m, MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyb, mlr_alloc_string_from_double(b, MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyn, mlr_alloc_string_from_ll(pstate->count), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyq, mlr_alloc_string_from_double(q, MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
if (pstate->do_verbose) {
lrec_put(poutrec, keyl1, mlr_alloc_string_from_double(l1, MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyl2, mlr_alloc_string_from_double(l2, MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyv11, mlr_alloc_string_from_double(v1[0], MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyv12, mlr_alloc_string_from_double(v1[1], MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyv21, mlr_alloc_string_from_double(v2[0], MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
lrec_put(poutrec, keyv22, mlr_alloc_string_from_double(v2[1], MLR_GLOBALS.ofmt), FREE_ENTRY_VALUE);
}
}
} else {
char* key = (pstate->do_which == DO_CORR) ? pstate->corr_output_field_name : pstate->cov_output_field_name;
if (pstate->count < 2LL) {
lrec_put(poutrec, key, "", NO_FREE);
} else {
double output = mlr_get_cov(pstate->count, pstate->sumx, pstate->sumy, pstate->sumxy);
if (pstate->do_which == DO_CORR) {
double sigmax = sqrt(mlr_get_var(pstate->count, pstate->sumx, pstate->sumx2));
double sigmay = sqrt(mlr_get_var(pstate->count, pstate->sumy, pstate->sumy2));
output = output / sigmax / sigmay;
}
char* val = mlr_alloc_string_from_double(output, MLR_GLOBALS.ofmt);
lrec_put(poutrec, key, val, FREE_ENTRY_VALUE);
}
}
}
static void linreg_pca_fit(void* pvstate, double x, double y, lrec_t* poutrec) {
stats2_corr_cov_state_t* pstate = pvstate;
if (!pstate->fit_ready) {
double Q[2][2];
mlr_get_cov_matrix(pstate->count,
pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy, Q);
double l1, l2; // Eigenvalues
double v1[2], v2[2]; // Eigenvectors
mlr_get_real_symmetric_eigensystem(Q, &l1, &l2, v1, v2);
double x_mean = pstate->sumx / pstate->count;
double y_mean = pstate->sumy / pstate->count;
mlr_get_linear_regression_pca(l1, l2, v1, v2, x_mean, y_mean, &pstate->m, &pstate->b, &pstate->q);
pstate->fit_ready = TRUE;
}
if (pstate->count < 2LL) {
lrec_put(poutrec, pstate->pca_fit_output_field_name, "", NO_FREE);
} else {
double yfit = pstate->m * x + pstate->b;
lrec_put(poutrec, pstate->pca_fit_output_field_name, mlr_alloc_string_from_double(yfit, MLR_GLOBALS.ofmt),
FREE_ENTRY_VALUE);
}
}
static void stats2_corr_cov_free(stats2_acc_t* pstats2_acc) {
stats2_corr_cov_state_t* pstate = pstats2_acc->pvstate;
free(pstate->covx_00_output_field_name);
free(pstate->covx_01_output_field_name);
free(pstate->covx_10_output_field_name);
free(pstate->covx_11_output_field_name);
free(pstate->pca_m_output_field_name);
free(pstate->pca_b_output_field_name);
free(pstate->pca_n_output_field_name);
free(pstate->pca_q_output_field_name);
free(pstate->pca_l1_output_field_name);
free(pstate->pca_l2_output_field_name);
free(pstate->pca_v11_output_field_name);
free(pstate->pca_v12_output_field_name);
free(pstate->pca_v21_output_field_name);
free(pstate->pca_v22_output_field_name);
free(pstate->pca_fit_output_field_name);
free(pstate->corr_output_field_name);
free(pstate->cov_output_field_name);
free(pstate);
free(pstats2_acc);
}
static stats2_acc_t* stats2_corr_cov_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name,
bivar_measure_t do_which, int do_verbose)
{
stats2_acc_t* pstats2_acc = mlr_malloc_or_die(sizeof(stats2_acc_t));
stats2_corr_cov_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_corr_cov_state_t));
pstate->count = 0LL;
pstate->sumx = 0.0;
pstate->sumy = 0.0;
pstate->sumx2 = 0.0;
pstate->sumxy = 0.0;
pstate->sumy2 = 0.0;
pstate->do_which = do_which;
pstate->do_verbose = do_verbose;
char* name1 = value_field_name_1;
char* name2 = value_field_name_2;
pstate->covx_00_output_field_name = mlr_paste_4_strings(name1, "_", name1, "_covx");
pstate->covx_01_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_covx");
pstate->covx_10_output_field_name = mlr_paste_4_strings(name2, "_", name1, "_covx");
pstate->covx_11_output_field_name = mlr_paste_4_strings(name2, "_", name2, "_covx");
pstate->pca_m_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_m");
pstate->pca_b_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_b");
pstate->pca_n_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_n");
pstate->pca_q_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_quality");
pstate->pca_l1_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_eival1");
pstate->pca_l2_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_eival2");
pstate->pca_v11_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_eivec11");
pstate->pca_v12_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_eivec12");
pstate->pca_v21_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_eivec21");
pstate->pca_v22_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_eivec22");
pstate->pca_fit_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_pca_fit");
pstate->fit_ready = FALSE;
pstate->m = -999.0;
pstate->b = -999.0;
pstate->corr_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_corr");
pstate->cov_output_field_name = mlr_paste_4_strings(name1, "_", name2, "_cov");
pstats2_acc->pvstate = (void*)pstate;
pstats2_acc->pingest_func = stats2_corr_cov_ingest;
pstats2_acc->pemit_func = stats2_corr_cov_emit;
if (do_which == DO_LINREG_PCA)
pstats2_acc->pfit_func = linreg_pca_fit;
else
pstats2_acc->pfit_func = NULL;
pstats2_acc->pfree_func = stats2_corr_cov_free;
return pstats2_acc;
}
static stats2_acc_t* stats2_corr_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
return stats2_corr_cov_alloc(value_field_name_1, value_field_name_2, stats2_acc_name, DO_CORR, do_verbose);
}
static stats2_acc_t* stats2_cov_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
return stats2_corr_cov_alloc(value_field_name_1, value_field_name_2, stats2_acc_name, DO_COV, do_verbose);
}
static stats2_acc_t* stats2_covx_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
return stats2_corr_cov_alloc(value_field_name_1, value_field_name_2, stats2_acc_name, DO_COVX, do_verbose);
}
static stats2_acc_t* stats2_linreg_pca_alloc(char* value_field_name_1, char* value_field_name_2, char* stats2_acc_name, int do_verbose) {
return stats2_corr_cov_alloc(value_field_name_1, value_field_name_2, stats2_acc_name, DO_LINREG_PCA, do_verbose);
}
Reference in a new issue View git blame Copy permalink