/*
 * linear_regression_smpl.c - linear_regression output sampling module for all PMU  models
 * Qiang Wu, Computer Architecture Group, Department of Computing, Imperial College London.
 * <qiangwu@doc.ic.ac.uk>
 *
 * Modified from perfmon_dfl_smpl.c
 *
 * Copyright (c) 2002-2006 Hewlett-Packard Development Company, L.P.
 * Contributed by Stephane Eranian <eranian@hpl.hp.com>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307 USA
 */
#include <math.h>
#include <stdlib.h>
//#include <asm/msr.h>  // for rdtscll

#include "pfmon.h"

#include <perfmon/perfmon_dfl_smpl.h>

#define SMPL_MOD_NAME			"linear-regression"

/*
 * forward declaration
 */
//pfmon_smpl_module_t linear_regression_smpl_module;

/*
 * types and defines
 */

/* Time stamps struct and functions to record time stamps elapsed */

#define rdtscll(val) do { \
unsigned int __a,__d; \
__asm__ __volatile__("rdtsc" : "=a" (__a), "=d" (__d)); \
(val) = ((unsigned long)__a) | (((unsigned long)__d)<<32); \
} while(0)

typedef struct timing_data_t {
	uint64_t total_ts;
	uint64_t smpl_ts;
} TIMING_DATA;

TIMING_DATA g_timestamps = { 0, 0 };

inline uint64_t getTotalTimeStamps( TIMING_DATA *pts )
{
	return pts->total_ts;
}

inline uint64_t getSmplTimeStamps( TIMING_DATA *pts )
{
	return pts->smpl_ts;
}

inline void initTimeStamps(struct timing_data_t *pts)
{
	pts->total_ts = 0;
	pts->smpl_ts = 0;
}

inline void startTotalTimeStamp(struct timing_data_t *pts, uint64_t start_ts)
{
	pts->total_ts = start_ts;
}

inline uint64_t endTotalTimeStamp(struct timing_data_t *pts, uint64_t end_ts)
{
	pts->total_ts = end_ts - pts->total_ts;

	return pts->total_ts;
}

inline void updateSmplTimeStamp(struct timing_data_t *pts, uint64_t ts)
{
	pts->smpl_ts += ts;
}

inline long double calcTimeOverhead( struct timing_data_t *pts )
{
	return (long double) 1.0 * pts->smpl_ts / pts->total_ts;
}


/* Options specific to linear regression module. */
typedef struct linear_regression_options_t {
	long double alpha;
	int printrecords;
	int recordselftime;
	unsigned long ycol;
} LINEAR_REGRESSION_OPTIONS;

LINEAR_REGRESSION_OPTIONS g_linear_regression_options = { 0.01, 0, 0, 0 };

#define LINEAR_REGRESSION_OPT_ALPHA_IDX	5
#define LINEAR_REGRESSION_OPT_YCOL_IDX	10

static struct option g_linear_regression_options_info[]={
	{ "linear-regression-alpha", 1, NULL, PFMON_OPT_SMPL_BASE + LINEAR_REGRESSION_OPT_ALPHA_IDX },
	{ "linear-regression-printrecords", 0, &(g_linear_regression_options.printrecords), 1 },
	{ "linear-regression-recordselftime", 0, &(g_linear_regression_options.recordselftime), 1 },
	{ "linear-regression-eventindex", 1, NULL, PFMON_OPT_SMPL_BASE + LINEAR_REGRESSION_OPT_YCOL_IDX },
	{ 0, 0, 0, 0}
};

inline long double getAlpha( LINEAR_REGRESSION_OPTIONS * popts )
{
	return popts->alpha;
}

inline void setAlpha( LINEAR_REGRESSION_OPTIONS * popts, long double a )
{
	popts->alpha = a;
}

inline unsigned long  getYCol( LINEAR_REGRESSION_OPTIONS * popts )
{
	return popts->ycol;
}

inline void setYCol( LINEAR_REGRESSION_OPTIONS * popts, unsigned long n )
{
	popts->ycol = n;
}

inline int doesPrintRecords( LINEAR_REGRESSION_OPTIONS * popts )
{
	return popts->printrecords;
}

inline int doesRecordSelfTime( LINEAR_REGRESSION_OPTIONS * popts )
{
	return popts->recordselftime;
}

inline void initLinearRegressionOptions(struct linear_regression_options_t *popts)
{
	popts->alpha = 0.01;
	popts->printrecords = 0;
	popts->recordselftime = 0;
	popts->ycol = 0;
}

static void linear_regression_show_options(void)
{
	printf(	"\t--linear-regression-alpha=value		Set the alpha to value (default is 0.01)\n"
		"\t--linear-regression-printrecords		Print the records sampled (default no)\n"
		"\t--linear-regression-recordselftime		Record and show the time consumed by the sampling module itself (default no)\n"
		"\t--linear-regression-eventindex=number	Event index after the sampling event to do linear regression (default 0, the one just after sampling event)\n"
		"\t                                 		If the given event index is too large, the modular against number of events excluding the sampling\n"
		"\t                              		event will be used.\n"
	);
}

static int linear_regression_parse_options(int code, char *optarg)
{
	switch(code) {
		case (PFMON_OPT_SMPL_BASE + LINEAR_REGRESSION_OPT_ALPHA_IDX):
			setAlpha( &g_linear_regression_options, (long double) strtold(optarg, NULL) );
			break;
		case (PFMON_OPT_SMPL_BASE + LINEAR_REGRESSION_OPT_YCOL_IDX):
			setYCol( &g_linear_regression_options, strtoul(optarg, NULL, 0) );
			break;
		default:
			return -1;
	}

	return 0;
}


/* Function to print sampling record. */

inline int printRecord(FILE *fp, uint64_t entry, pfm_dfl_smpl_entry_t *ent, uint64_t *reg, pfmon_event_set_t *active_set, unsigned int npmds)
{
	unsigned int n;
	int ret;

	ret = fprintf(	fp, "%-8"PRIu64" %8d %8d %2d 0x%llx 0x%016"PRIx64" %3u %5u %16"PRId64" ",
			entry,
			ent->pid,
			ent->tgid,
			ent->cpu,
			(unsigned long long) ent->ip,
			ent->tstamp,
			(unsigned int) ent->ovfl_pmd,
			ent->set,
			ent->last_reset_val);

	for (n = 0; n < npmds; n++) {
		ret = fprintf(fp, " %"PRIu64"", reg[active_set->rev_smpl_pmds[ent->ovfl_pmd].map_pmd_evt[n].off]);
	}

	ret = fputc('\n', fp);

	return ret;
}

/* Line parameters struct and functions to init, start, update and print it. */
typedef struct line_params_t {
	unsigned long num_lines;
	unsigned long num_points;
	long double start_x;
	long double start_y;
	long double end_x;
	long double end_y;
	long double slope;
	long double intercept;
	long double sum_x; 
	long double sum_y;
	long double sum_xx;
	long double sum_yy;
	long double sum_xy;
	long double residuals;
	long double sigma;
	long double max_sigma;
	long double init_y;
} LINE_PARAMETERS;

LINE_PARAMETERS g_line = { 
	0,
	0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0,
	0.0
};

inline unsigned long getNumLines( LINE_PARAMETERS * pline )
{
	return pline->num_lines;
}

inline unsigned long getNumPoints( LINE_PARAMETERS * pline )
{
	return pline->num_points;
}

inline unsigned long getLineEndX( LINE_PARAMETERS * pline )
{
	return pline->end_x;
}

inline long double getLineSigma( LINE_PARAMETERS * pline )
{
	return pline->sigma;
}

inline void initLine(LINE_PARAMETERS *pline)
{
	pline->num_lines = 0;
	pline->num_points = 0;
	pline->start_x  = 0.0;
	pline->start_y  = 0.0;
	pline->end_x  = 0.0;
	pline->end_y  = 0.0;
	pline->slope = 0.0;
	pline->intercept = 0.0;
	pline->sum_x = 0.0;
	pline->sum_y = 0.0;
	pline->sum_xx = 0.0;
	pline->sum_yy = 0.0;
	pline->sum_xy = 0.0;
	pline->residuals = 0.0;
	pline->sigma = 0.0;
	pline->max_sigma = 0.0;
}

void startLine(LINE_PARAMETERS *pline, long double x, long double y)
{
	if ( pline->num_points == 0 ) {
		pline->num_points = 1;
		pline->start_x = x;
		pline->start_y = y;
		pline->end_x = x;
		pline->end_y = y;
		pline->sum_x = x;
		pline->sum_y = y;
		pline->sum_xx = x * x;
		pline->sum_yy = y * y;
		pline->sum_xy = x * y;
		pline->slope = 0.0;
		pline->intercept = 0.0;
		pline->residuals = 0.0;
		pline->sigma = 0.0;
	}
	else {
		pline->num_points = 2;
		pline->start_x = pline->end_x;
		pline->start_y = pline->end_y;
		pline->end_x = x;
		pline->end_y = y;
		pline->sum_x = pline->start_x + x;
		pline->sum_y = pline->start_y + y;
		pline->sum_xx = pline->start_x * pline->start_x + x * x;
		pline->sum_yy = pline->start_y * pline->start_y + y * y;
		pline->sum_xy = pline->start_x * pline->start_y + x * y;
		// FIXME: risk of divide 0 if start_x == end_x
		pline->slope = ( pline->end_y - pline->start_y ) / ( pline->end_x - pline->start_x );
		pline->intercept = ( pline->start_y * pline->end_x - pline->start_x * pline->end_y ) / ( pline->end_x - pline->start_x );
		pline->residuals = 0.0;
		pline->sigma = 0.0;
	}

	if (pline->num_lines == 0) { pline->init_y = y; }

	pline->num_lines ++;
}

void updateLine(LINE_PARAMETERS *pline, long double new_x, long double new_y)
{
	pline->num_points ++;
	if ( pline->num_points < 2 ) {
		pline->start_x = new_x;
		pline->start_y = new_y;
		pline->end_x = new_x;
		pline->end_y = new_y;
		pline->sum_x = new_x;
		pline->sum_y = new_y;
		pline->sum_xx = new_x * new_x;
		pline->sum_yy = new_y * new_y;
		pline->sum_xy = new_x * new_y;
		pline->slope = 0.0;
		pline->intercept = 0.0;
		pline->residuals = 0.0;
		pline->sigma = 0.0;
	}
	else {
		pline->end_x = new_x;
		pline->end_y = new_y;
		pline->sum_x += new_x;
		pline->sum_y += new_y;
		pline->sum_xx += new_x * new_x;
		pline->sum_yy += new_y * new_y;
		pline->sum_xy += new_x * new_y;
		pline->slope = ( pline->num_points * pline->sum_xy - pline->sum_x * pline->sum_y ) / ( pline->num_points * pline->sum_xx - pline->sum_x * pline->sum_x );
		pline->intercept = ( pline->sum_xx * pline->sum_y - pline->sum_xy * pline->sum_x ) / ( pline->num_points * pline->sum_xx - pline->sum_x * pline->sum_x );
		pline->residuals = (pline->sum_yy - 2 * pline->intercept * pline->sum_y + pline->slope * pline->slope * pline->sum_xx \
				- 2 * pline->slope * pline->sum_xy + pline->num_points * pline->intercept * pline->intercept \
				+ 2 * pline->slope * pline->intercept * pline->sum_x);
		pline->sigma = sqrt( pline->residuals / (pline->num_points -2) );
		if (pline->sigma > pline->max_sigma) { pline->max_sigma = pline->sigma; }
	}
}

inline int out3SigmaFromLine(LINE_PARAMETERS *pline, long double alpha, long double x, long double y)
{
	long double y_hat;

	y_hat = pline->slope * x + pline->intercept;

	if (pline->sigma > 0.0) {
		return fabs(y - y_hat) > 3 * pline->sigma;
	}
	else {
		return fabs(y - y_hat) > alpha * fabs(y_hat);
	}
}

inline long double calcNormalizedMaxSigma( LINE_PARAMETERS * pline)
{
	return pline->max_sigma / fabs(pline->end_y - pline->init_y);
}

inline int printLine(FILE *fp, LINE_PARAMETERS *pline)
{
	// LINE=(start_x, start_y, end_x, end_y, slope, intercept, num_points, sigma)
	
	return fprintf( fp, "LINE#%lu = (%lu, %Lf, %Lf, %Lf, %Lf, %Lf, %Lf, %Lf)\n", pline->num_lines, pline->num_points,
			pline->start_x, pline->start_y, pline->end_x, pline->end_y, pline->slope, pline->intercept, pline->sigma);
}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*
 * functions of sample module
 */
static int
linear_regression_process_samples(pfmon_sdesc_t *sdesc)
{
	pfm_dfl_smpl_hdr_t *hdr;
	pfm_dfl_smpl_entry_t *ent;
	pfmon_smpl_desc_t *csmpl;
	pfmon_event_set_t *active_set;
	uint64_t skip, count, entry;
	FILE *fp;
	uint64_t *reg;
	unsigned int npmds, n;
	int ret;

	long double x, y;

	unsigned long start_ts, end_ts;

	rdtscll(start_ts);

 	csmpl = &sdesc->csmpl;
	hdr   = csmpl->smpl_hdr;
 	fp = csmpl->smpl_fp;

	if (hdr == NULL)
		return -1;

	ent	   = (pfm_dfl_smpl_entry_t *)(hdr+1);
	entry      = options.opt_aggr ? *csmpl->aggr_count : csmpl->entry_count;
	count      = hdr->hdr_count;
	active_set = sdesc->sets; /* only one set when sampling */

	DPRINT(("count=%"PRIu64" entry=%"PRIu64"\n", count, entry));

	/*
 	 * check if we have new entries
 	 * if so skip the old entries and process only the new ones
 	 */
	if((csmpl->last_ovfl == hdr->hdr_overflows && csmpl->last_count <= count)
	  || ((csmpl->last_ovfl+1) == hdr->hdr_overflows && csmpl->last_count < count)) {
		skip = csmpl->last_count;
		vbprintf("skip %"PRIu64" samples out of %"PRIu64" (overflows: %"PRIu64")\n",
			  skip,
			  count,
			  hdr->hdr_overflows);
	} else {
		skip = 0;
	}

	vbprintf("skip=%"PRIu64" count=%"PRIu64" ovfl=%"PRIu64" last_o=%"PRIu64" last_c=%"PRIu64"\n",
		skip,
		count,
		hdr->hdr_overflows,
		csmpl->last_ovfl,
		csmpl->last_count);

	/*
 	 * only account for new entries, i.e., skip leftover entries
 	 */
	if (options.opt_aggr) {
		*csmpl->aggr_count += count - skip;
	} else {
		csmpl->entry_count += count - skip;
	}

	csmpl->last_count = count;
	csmpl->last_ovfl = hdr->hdr_overflows;

	while(count--) {
		reg = (uint64_t *)(ent+1);

		/* XXX: why not put npmds outside the while loop? not clear about pfmon implementation idea... */
		npmds = active_set->rev_smpl_pmds[ent->ovfl_pmd].num_smpl_pmds;

		if (skip == 0) {
			if ( doesPrintRecords(&g_linear_regression_options) ) {
				printRecord( fp, entry, ent, reg, active_set, npmds );
			}

			x = getLineEndX(&g_line) + (-1.0) * ((int64_t) ent->last_reset_val);
			if ( npmds > 0 ) {
				n = getYCol(&g_linear_regression_options) % npmds;
				y = reg[active_set->rev_smpl_pmds[ent->ovfl_pmd].map_pmd_evt[n].off];
			}
			else {
				/* use the time stamp value as y. */
				y = ent->tstamp;
			}

			/* check if at the very beginning */
			if ( getNumPoints(&g_line) == 0 ) {
				startLine(&g_line, x, y);
			}
			else if ( getNumPoints(&g_line) > 1 && out3SigmaFromLine( &g_line, getAlpha(&g_linear_regression_options), x, y ) ) {
				ret = printLine( fp, &g_line );
				if (ret == -1) { goto error; }
				startLine(&g_line, x, y);
			}
			else {
				updateLine(&g_line, x, y);
			}

			entry++;
		}
		else {
			skip--;
		}

		reg += npmds;

		ent  = (pfm_dfl_smpl_entry_t *)reg;
	}

	if ( doesRecordSelfTime(&g_linear_regression_options) ) {
		rdtscll(end_ts);
		updateSmplTimeStamp( &g_timestamps, end_ts - start_ts );
	}

	return 0;
error:
	warning("cannot write to sampling file: %s\n", strerror(errno));
	if (doesRecordSelfTime(&g_linear_regression_options)) {
		rdtscll(end_ts);
		updateSmplTimeStamp( &g_timestamps, end_ts - start_ts );
	}
	return -1;
}


static int
linear_regression_print_header(pfmon_sdesc_t *sdesc)
{
	pfmon_smpl_desc_t *csmpl = &sdesc->csmpl;
	FILE *fp = csmpl->smpl_fp;

	fprintf( fp, "# linear regression options: Alpha=%Lf, PrintRecord=%d, SelfTiming=%d, YCol=%lu\n#\n", 
			getAlpha(&g_linear_regression_options),
			doesPrintRecords(&g_linear_regression_options),
			doesRecordSelfTime(&g_linear_regression_options),
			getYCol(&g_linear_regression_options));

	fprintf(fp, "# description of LINE output:\n"
		    "#\tLINE#<line_id>=(<number_of_points>, <start_x>, <start_y>, <end_x>, <end_y>, <slope>, <intercept>, <estimated_standard_deviation>\n#\n"
		);

	if ( doesPrintRecords(&g_linear_regression_options) ) {
		fprintf(fp, "# description of RECORD columns:\n"
			    "#\tcolumn  1: entry number\n"
		 	    "#\tcolumn  2: process id\n"
		 	    "#\tcolumn  3: thread id\n"
			    "#\tcolumn  4: cpu number\n"
			    "#\tcolumn  5: instruction pointer\n"
			    "#\tcolumn  6: unique timestamp\n"
			    "#\tcolumn  7: overflowed PMD index\n"
			    "#\tcolumn  8: event set\n"
			    "#\tcolumn  9: initial value of overflowed PMD (sampling period)\n"
			    "#\tfollowed by optional sampled PMD values in command line order\n");
	}

	return 0;
}

static int linear_regression_initialize_module(void)
{
	initLinearRegressionOptions( &g_linear_regression_options );

	return pfmon_register_options(g_linear_regression_options_info, sizeof(g_linear_regression_options_info));
}


pfmon_smpl_module_t linear_regression_smpl_module;

static void linear_regression_initialize_mask(void)
{
	pfmon_bitmask_setall(linear_regression_smpl_module.pmu_mask);
}

static int
linear_regression_initialize_session(pfmon_sdesc_t *sdesc)
{
	uint64_t ts;

	if (doesRecordSelfTime(&g_linear_regression_options)) {
		initTimeStamps( &g_timestamps );

		rdtscll(ts);
		startTotalTimeStamp( &g_timestamps, ts );
	}

	initLine( &g_line );

	return 0;
}

static int
linear_regression_terminate_session(pfmon_sdesc_t *sdesc)
{
	pfmon_smpl_desc_t *csmpl;

	uint64_t end_ts;
	uint64_t total_ts;
	long double overhead;
	long double nmesd;

	csmpl = &sdesc->csmpl;

	nmesd = calcNormalizedMaxSigma( &g_line );

	printLine( csmpl->smpl_fp, &g_line );

	if (doesRecordSelfTime(&g_linear_regression_options)) {
		rdtscll(end_ts);
		total_ts = endTotalTimeStamp( &g_timestamps, end_ts );
		overhead = calcTimeOverhead( &g_timestamps );

		fprintf(csmpl->smpl_fp, "# total time stamps (from initialize_session to terminate_session) : %"PRIu64"\n", total_ts);
		fprintf(csmpl->smpl_fp, "# sampling time stamps : %lu\n", getSmplTimeStamps(&g_timestamps));
		fprintf(csmpl->smpl_fp, "# time overhead : %Lf\n", overhead);
	}

	fprintf(csmpl->smpl_fp, "# total number of samples : %"PRIu64"\n", csmpl->entry_count);
	fprintf(csmpl->smpl_fp, "# total number of lines : %lu\n", getNumLines(&g_line));
	fprintf(csmpl->smpl_fp, "# ratio of samples to lines : %f\n", 1.0 * csmpl->entry_count / getNumLines(&g_line));
	fprintf(csmpl->smpl_fp, "# normalized maximal estimated standard deviation among lines : %Lf\n", nmesd);
	//fprintf(csmpl->smpl_fp, "# number of buffer overflows : %"PRIu64"\n", csmpl->last_ovfl);

	return 0;
}

pfmon_smpl_module_t linear_regression_smpl_module = {
	.name		    = SMPL_MOD_NAME,
	.description	    = "Linear regression of sample values",
	.initialize_module  = linear_regression_initialize_module,
	.initialize_mask    = linear_regression_initialize_mask,
	.initialize_session = linear_regression_initialize_session,
	.terminate_session  = linear_regression_terminate_session,
	.process_samples    = linear_regression_process_samples,
	.print_header       = linear_regression_print_header,
	.parse_options      = linear_regression_parse_options,
	.show_options       = linear_regression_show_options,
	.fmt_name	    = PFM_DFL_SMPL_NAME
};