<commit_before>#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "lib/mlrutil.h"
#include "lib/mlrstat.h"
#include "containers/sllv.h"
#include "containers/slls.h"
#include "containers/lhmslv.h"
#include "containers/lhms2v.h"
#include "containers/lhmsv.h"
#include "containers/mixutil.h"
#include "mapping/mappers.h"
#include "cli/argparse.h"

#define DO_CORR 0x11
#define DO_COV  0x22
#define DO_COVX 0x33

// ----------------------------------------------------------------
typedef void stats2_put_func_t(void* pvstate, double x, double y);
typedef void stats2_get_func_t(void* pvstate, char* name1, char* name2, lrec_t* poutrec);

typedef struct _stats2_t {
	void* pvstate;
	stats2_put_func_t* pput_func;
	stats2_get_func_t* pget_func;
} stats2_t;

typedef stats2_t* stats2_alloc_func_t(static_context_t* pstatx);

// ----------------------------------------------------------------
// Univariate linear regression
// ----------------------------------------------------------------
// There are N (xi, yi) pairs.
//
// E = sum (yi - m xi - b)^2
//
// DE/Dm = sum 2 (yi - m xi - b) (-xi) = 0
// DE/Db = sum 2 (yi - m xi - b) (-1)  = 0
//
// sum (yi - m xi - b) (xi) = 0
// sum (yi - m xi - b)      = 0
//
// sum (xi yi - m xi^2 - b xi) = 0
// sum (yi - m xi - b)         = 0
//
// m sum(xi^2) + b sum(xi) = sum(xi yi)
// m sum(xi)   + b N       = sum(yi)
//
// [ sum(xi^2)   sum(xi) ] [ m ] = [ sum(xi yi) ]
// [ sum(xi)     N       ] [ b ] = [ sum(yi)    ]
//
// [ m ] = [ sum(xi^2) sum(xi) ]^-1  [ sum(xi yi) ]
// [ b ]   [ sum(xi)   N       ]     [ sum(yi)    ]
//
//       = [ N         -sum(xi)  ]  [ sum(xi yi) ] * 1/D
//         [ -sum(xi)   sum(xi^2)]  [ sum(yi)    ]
//
// where
//
//   D = N sum(xi^2) - sum(xi)^2.
//
// So
//
//      N sum(xi yi) - sum(xi) sum(yi)
// m = --------------------------------
//                   D
//
//      -sum(xi)sum(xi yi) + sum(xi^2) sum(yi)
// b = ----------------------------------------
//                   D

typedef struct _stats2_linreg_state_t {
	unsigned long long count;
	double sumx;
	double sumy;
	double sumx2;
	double sumxy;
	static_context_t* pstatx;
} stats2_linreg_state_t;
void stats2_linreg_put(void* pvstate, double x, double y) {
	stats2_linreg_state_t* pstate = pvstate;
	pstate->count++;
	pstate->sumx  += x;
	pstate->sumy  += y;
	pstate->sumx2 += x*x;
	pstate->sumxy += x*y;
}
void stats2_linreg_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
	stats2_linreg_state_t* pstate = pvstate;
	int n = pstate->count;
	double sumx  = pstate->sumx;
	double sumy  = pstate->sumy;
	double sumx2 = pstate->sumx2;
	double sumxy = pstate->sumxy;
	double D =  n * sumx2 - sumx*sumx;
	double m = (n * sumxy - sumx * sumy) / D;
	double b = (-sumx * sumxy + sumx2 * sumy) / D;
	// xxx gah ... need a 2nd pass through the data to get the error-bars.

	// xxx make a 2nd filter to compute the error-bars given the data & the m & the b?
	//
	//	# Young 1962, pp. 122-124.  Compute sample variance of linear
	//	# approximations, then variances of m and b.
	//	var_z = 0.0
	//	for i in range(0, N):
	//		var_z += (m * xs[i] + b - ys[i])**2
	//	var_z /= N
	//
	//	var_m = (N * var_z) / D
	//	var_b = (var_z * sumx2) / D
	//
	//	return [m, b, math.sqrt(var_m), math.sqrt(var_b)]

	char* key = mlr_paste_4_strings(name1, "_", name2, "_m");
	char* val = mlr_alloc_string_from_double(m, pstate->pstatx->ofmt);
	lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);

	key = mlr_paste_4_strings(name1, "_", name2, "_b");
	val = mlr_alloc_string_from_double(b, pstate->pstatx->ofmt);
	lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
}
stats2_t* stats2_linreg_alloc(static_context_t* pstatx) {
	stats2_t* pstats2 = mlr_malloc_or_die(sizeof(stats2_t));
	stats2_linreg_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_linreg_state_t));
	pstate->count = 0LL;
	pstate->sumx  = 0.0;
	pstate->sumy  = 0.0;
	pstate->sumx2 = 0.0;
	pstate->sumxy = 0.0;
	pstate->pstatx = pstatx;
	pstats2->pvstate = (void*)pstate;
	pstats2->pput_func = &stats2_linreg_put;
	pstats2->pget_func = &stats2_linreg_get;
	return pstats2;
}

// ----------------------------------------------------------------
// 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;
	static_context_t* pstatx;
} stats2_r2_state_t;
void stats2_r2_put(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;
}
void stats2_r2_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
	stats2_r2_state_t* pstate = pvstate;
	char* suffix = "r2";
	char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix);
	if (pstate->count < 2LL) {
		lrec_put(poutrec, key, "", LREC_FREE_ENTRY_KEY);
	} 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, pstate->pstatx->ofmt);
		lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
	}
}
stats2_t* stats2_r2_alloc(static_context_t* pstatx) {
	stats2_t* pstats2 = mlr_malloc_or_die(sizeof(stats2_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->pstatx = pstatx;
	pstats2->pvstate   = (void*)pstate;
	pstats2->pput_func = &stats2_r2_put;
	pstats2->pget_func = &stats2_r2_get;
	return pstats2;
}

// ----------------------------------------------------------------
// def find_sample_covariance(xs, ys):
//      N = len(xs)
//      mean_x = find_mean(xs)
//      mean_y = find_mean(ys)
//
//      sum = 0.0
//      for k in range(0, N):
//              sum += (xs[k] - mean_x) * (ys[k] - mean_y)
//
//      return sum / (N-1.0)

// 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;
	int    do_which;
	static_context_t* pstatx;
} stats2_corr_cov_state_t;
void stats2_corr_cov_put(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;
}

void stats2_corr_cov_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
	stats2_corr_cov_state_t* pstate = pvstate;
	if (pstate->do_which == DO_COVX) {
		char* key00 = mlr_paste_4_strings(name1, "_", name2, "_q00");
		char* key01 = mlr_paste_4_strings(name1, "_", name2, "_q01");
		char* key10 = mlr_paste_4_strings(name1, "_", name2, "_q10");
		char* key11 = mlr_paste_4_strings(name1, "_", name2, "_q11");
		if (pstate->count < 2LL) {
			lrec_put(poutrec, key00, "", LREC_FREE_ENTRY_KEY);
			lrec_put(poutrec, key01, "", LREC_FREE_ENTRY_KEY);
			lrec_put(poutrec, key10, "", LREC_FREE_ENTRY_KEY);
			lrec_put(poutrec, key11, "", LREC_FREE_ENTRY_KEY);
		} else {
			double q00, q01, q10, q11;
			mlr_get_cov_matrix(pstate->count,
				pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy,
				&q00, &q01, &q10, &q11);
			char* val00 = mlr_alloc_string_from_double(q00, pstate->pstatx->ofmt);
			char* val01 = mlr_alloc_string_from_double(q01, pstate->pstatx->ofmt);
			char* val10 = mlr_alloc_string_from_double(q10, pstate->pstatx->ofmt);
			char* val11 = mlr_alloc_string_from_double(q11, pstate->pstatx->ofmt);
			lrec_put(poutrec, key00, val00, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
			lrec_put(poutrec, key01, val01, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
			lrec_put(poutrec, key10, val10, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
			lrec_put(poutrec, key11, val11, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
		}
	} else {
		char* suffix = (pstate->do_which == DO_CORR) ? "corr" : "cov";
		char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix);
		if (pstate->count < 2LL) {
			lrec_put(poutrec, key, "", LREC_FREE_ENTRY_KEY);<commit_msg>nicer names for covariance-matrix elements<commit_after>#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "lib/mlrutil.h"
#include "lib/mlrstat.h"
#include "containers/sllv.h"
#include "containers/slls.h"
#include "containers/lhmslv.h"
#include "containers/lhms2v.h"
#include "containers/lhmsv.h"
#include "containers/mixutil.h"
#include "mapping/mappers.h"
#include "cli/argparse.h"

#define DO_CORR 0x11
#define DO_COV  0x22
#define DO_COVX 0x33

// ----------------------------------------------------------------
typedef void stats2_put_func_t(void* pvstate, double x, double y);
typedef void stats2_get_func_t(void* pvstate, char* name1, char* name2, lrec_t* poutrec);

typedef struct _stats2_t {
	void* pvstate;
	stats2_put_func_t* pput_func;
	stats2_get_func_t* pget_func;
} stats2_t;

typedef stats2_t* stats2_alloc_func_t(static_context_t* pstatx);

// ----------------------------------------------------------------
// Univariate linear regression
// ----------------------------------------------------------------
// There are N (xi, yi) pairs.
//
// E = sum (yi - m xi - b)^2
//
// DE/Dm = sum 2 (yi - m xi - b) (-xi) = 0
// DE/Db = sum 2 (yi - m xi - b) (-1)  = 0
//
// sum (yi - m xi - b) (xi) = 0
// sum (yi - m xi - b)      = 0
//
// sum (xi yi - m xi^2 - b xi) = 0
// sum (yi - m xi - b)         = 0
//
// m sum(xi^2) + b sum(xi) = sum(xi yi)
// m sum(xi)   + b N       = sum(yi)
//
// [ sum(xi^2)   sum(xi) ] [ m ] = [ sum(xi yi) ]
// [ sum(xi)     N       ] [ b ] = [ sum(yi)    ]
//
// [ m ] = [ sum(xi^2) sum(xi) ]^-1  [ sum(xi yi) ]
// [ b ]   [ sum(xi)   N       ]     [ sum(yi)    ]
//
//       = [ N         -sum(xi)  ]  [ sum(xi yi) ] * 1/D
//         [ -sum(xi)   sum(xi^2)]  [ sum(yi)    ]
//
// where
//
//   D = N sum(xi^2) - sum(xi)^2.
//
// So
//
//      N sum(xi yi) - sum(xi) sum(yi)
// m = --------------------------------
//                   D
//
//      -sum(xi)sum(xi yi) + sum(xi^2) sum(yi)
// b = ----------------------------------------
//                   D

typedef struct _stats2_linreg_state_t {
	unsigned long long count;
	double sumx;
	double sumy;
	double sumx2;
	double sumxy;
	static_context_t* pstatx;
} stats2_linreg_state_t;
void stats2_linreg_put(void* pvstate, double x, double y) {
	stats2_linreg_state_t* pstate = pvstate;
	pstate->count++;
	pstate->sumx  += x;
	pstate->sumy  += y;
	pstate->sumx2 += x*x;
	pstate->sumxy += x*y;
}
void stats2_linreg_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
	stats2_linreg_state_t* pstate = pvstate;
	int n = pstate->count;
	double sumx  = pstate->sumx;
	double sumy  = pstate->sumy;
	double sumx2 = pstate->sumx2;
	double sumxy = pstate->sumxy;
	double D =  n * sumx2 - sumx*sumx;
	double m = (n * sumxy - sumx * sumy) / D;
	double b = (-sumx * sumxy + sumx2 * sumy) / D;
	// xxx gah ... need a 2nd pass through the data to get the error-bars.

	// xxx make a 2nd filter to compute the error-bars given the data & the m & the b?
	//
	//	# Young 1962, pp. 122-124.  Compute sample variance of linear
	//	# approximations, then variances of m and b.
	//	var_z = 0.0
	//	for i in range(0, N):
	//		var_z += (m * xs[i] + b - ys[i])**2
	//	var_z /= N
	//
	//	var_m = (N * var_z) / D
	//	var_b = (var_z * sumx2) / D
	//
	//	return [m, b, math.sqrt(var_m), math.sqrt(var_b)]

	char* key = mlr_paste_4_strings(name1, "_", name2, "_m");
	char* val = mlr_alloc_string_from_double(m, pstate->pstatx->ofmt);
	lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);

	key = mlr_paste_4_strings(name1, "_", name2, "_b");
	val = mlr_alloc_string_from_double(b, pstate->pstatx->ofmt);
	lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
}
stats2_t* stats2_linreg_alloc(static_context_t* pstatx) {
	stats2_t* pstats2 = mlr_malloc_or_die(sizeof(stats2_t));
	stats2_linreg_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_linreg_state_t));
	pstate->count = 0LL;
	pstate->sumx  = 0.0;
	pstate->sumy  = 0.0;
	pstate->sumx2 = 0.0;
	pstate->sumxy = 0.0;
	pstate->pstatx = pstatx;
	pstats2->pvstate = (void*)pstate;
	pstats2->pput_func = &stats2_linreg_put;
	pstats2->pget_func = &stats2_linreg_get;
	return pstats2;
}

// ----------------------------------------------------------------
// 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;
	static_context_t* pstatx;
} stats2_r2_state_t;
void stats2_r2_put(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;
}
void stats2_r2_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
	stats2_r2_state_t* pstate = pvstate;
	char* suffix = "r2";
	char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix);
	if (pstate->count < 2LL) {
		lrec_put(poutrec, key, "", LREC_FREE_ENTRY_KEY);
	} 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, pstate->pstatx->ofmt);
		lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
	}
}
stats2_t* stats2_r2_alloc(static_context_t* pstatx) {
	stats2_t* pstats2 = mlr_malloc_or_die(sizeof(stats2_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->pstatx = pstatx;
	pstats2->pvstate   = (void*)pstate;
	pstats2->pput_func = &stats2_r2_put;
	pstats2->pget_func = &stats2_r2_get;
	return pstats2;
}

// ----------------------------------------------------------------
// def find_sample_covariance(xs, ys):
//      N = len(xs)
//      mean_x = find_mean(xs)
//      mean_y = find_mean(ys)
//
//      sum = 0.0
//      for k in range(0, N):
//              sum += (xs[k] - mean_x) * (ys[k] - mean_y)
//
//      return sum / (N-1.0)

// 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;
	int    do_which;
	static_context_t* pstatx;
} stats2_corr_cov_state_t;
void stats2_corr_cov_put(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;
}

void stats2_corr_cov_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
	stats2_corr_cov_state_t* pstate = pvstate;
	if (pstate->do_which == DO_COVX) {
		char* key00 = mlr_paste_4_strings(name1, "_", name1, "_covx");
		char* key01 = mlr_paste_4_strings(name1, "_", name2, "_covx");
		char* key10 = mlr_paste_4_strings(name2, "_", name1, "_covx");
		char* key11 = mlr_paste_4_strings(name2, "_", name2, "_covx");
		if (pstate->count < 2LL) {
			lrec_put(poutrec, key00, "", LREC_FREE_ENTRY_KEY);
			lrec_put(poutrec, key01, "", LREC_FREE_ENTRY_KEY);
			lrec_put(poutrec, key10, "", LREC_FREE_ENTRY_KEY);
			lrec_put(poutrec, key11, "", LREC_FREE_ENTRY_KEY);
		} else {
			double q00, q01, q10, q11;
			mlr_get_cov_matrix(pstate->count,
				pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy,
				&q00, &q01, &q10, &q11);
			char* val00 = mlr_alloc_string_from_double(q00, pstate->pstatx->ofmt);
			char* val01 = mlr_alloc_string_from_double(q01, pstate->pstatx->ofmt);
			char* val10 = mlr_alloc_string_from_double(q10, pstate->pstatx->ofmt);
			char* val11 = mlr_alloc_string_from_double(q11, pstate->pstatx->ofmt);
			lrec_put(poutrec, key00, val00, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
			lrec_put(poutrec, key01, val01, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
			lrec_put(poutrec, key10, val10, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
			lrec_put(poutrec, key11, val11, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
		}
	} else {
		char* suffix = (pstate->do_which == DO_CORR) ? "corr" : "cov";
		char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix);
		if (pstate->count < 2LL) {
			lrec_put(poutrec, key, "", LREC_FREE_ENTRY_KEY);
{'fraction_non_alphanumeric': 0.19645853827903478, 'fraction_numerical': 0.030553787396562698, 'mean_word_length': 2.9099547511312216, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 155, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5905755', 'n_tokens_mistral': 7768, 'n_tokens_neox': 7365, 'n_words': 1467}