improving multiple regression models

include/starpu_perfmodel.h

@@ -145,7 +145,8 @@ struct starpu_perfmodel
 	unsigned benchmarking;
 	unsigned is_init;
 
-	double *parameters;
+	void (*parameters)(struct starpu_task * task, double *parameters);
+	const char **parameters_names;
 	unsigned nparameters;
 	unsigned **combinations;
 	unsigned ncombinations;

src/core/perfmodel/multiple_regression.c

@@ -0,0 +1,264 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010, 2011, 2015  Université de Bordeaux
+ * Copyright (C) 2010, 2011  CNRS
+ *
+ * StarPU is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation; either version 2.1 of the License, or (at
+ * your option) any later version.
+ *
+ * StarPU 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 Lesser General Public License in COPYING.LGPL for more details.
+ */
+
+/* Code for computing multiple linear regression */
+
+#include <core/perfmodel/multiple_regression.h>
+
+typedef struct { int h, w; double *x;} matrix_t, *matrix;
+
+static void dump_multiple_regression_list(double *mx, double *my, unsigned ncoeff, unsigned nparameters, unsigned **combinations, struct starpu_perfmodel_history_list *list_history)
+{
+	struct starpu_perfmodel_history_list *ptr = list_history;
+	unsigned i = 0;
+
+	while (ptr)
+	{
+		mx[i*ncoeff] = 1.;
+		for(int j=0; j<ncoeff-1; j++)
+		{
+			mx[i*ncoeff+j+1] = 1.;
+			for(int k=0; k < nparameters; k++)
+				mx[i*ncoeff+j+1] *= pow(ptr->entry->parameters[k],combinations[j][k]);
+		}
+		my[i] = ptr->entry->duration;
+
+		ptr = ptr->next;
+		i++;
+	}
+}
+
+double dot(double *a, double *b, int len, int step)
+{
+	double r = 0;
+	while (len--) {
+		r += *a++ * *b;
+		b += step;
+	}
+	return r;
+}
+
+matrix mat_new(int h, int w)
+{
+	matrix r = malloc(sizeof(matrix_t) + sizeof(double) * w * h);
+	r->h = h, r->w = w;
+	r->x = (double*)(r + 1);
+	return r;
+}
+
+void mat_free(matrix a)
+{
+	free(a->x);
+	free(a);
+}
+
+matrix mat_mul(matrix a, matrix b)
+{
+	matrix r;
+	double *p, *pa;
+	int i, j;
+	if (a->w != b->h) return 0;
+
+	r = mat_new(a->h, b->w);
+	p = r->x;
+	for (pa = a->x, i = 0; i < a->h; i++, pa += a->w)
+		for (j = 0; j < b->w; j++)
+			*p++ = dot(pa, b->x + j, a->w, b->w);
+	return r;
+}
+
+void mat_show(matrix a)
+{
+	int i, j;
+	double *p = a->x;
+	for (i = 0; i < a->h; i++, putchar('\n'))
+		for (j = 0; j < a->w; j++)
+			printf("\t%7.3f", *p++);
+	putchar('\n');
+}
+
+// Inspired from: https://rosettacode.org/wiki/Matrix_transposition#C
+matrix transpose(matrix src)
+{
+	int i, j;
+	matrix dst;
+	dst = mat_new(src->w, src->h);
+
+	for (i = 0; i < src->h; i++)
+	  for (j = 0; j < src->w; j++)
+		dst->x[j * dst->w + i] = src->x[i * src->w + j];
+
+	return dst;
+}
+
+// Inspired from: http://www.programming-techniques.com/2011/09/numerical-methods-inverse-of-nxn-matrix.html
+matrix mat_inv(matrix src)
+{
+	int n = src->h;
+	int n2=2*n;
+    int i,j, k;
+	double ratio, a;
+	matrix r, dst;
+	r = mat_new(n, n2);
+	dst = mat_new(n, n);
+
+	for (i = 0; i < n; i++)
+	  for (j = 0; j < n; j++)
+		r->x[i*n2+j] = src->x[i*n+j];
+
+	for(i = 0; i < n; i++){
+	  for(j = n; j < 2*n; j++){
+            if(i==(j-n))
+                r->x[i*n2+j] = 1.0;
+            else
+                r->x[i*n2+j] = 0.0;
+	  }
+	}
+
+	for(i = 0; i < n; i++){
+	  for(j = 0; j < n; j++){
+	    if(i!=j){
+               for(k = 0; k < 2*n; k++){
+                   r->x[j*n2+k] -= (r->x[j*n2+i] / r->x[i*n2+i]) * r->x[i*n2+k];
+               }
+        }
+	  }
+	}
+
+	for(i = 0; i < n; i++){
+	  a = r->x[i*n2+i];
+	  for(j = 0; j < 2*n; j++){
+            r->x[i*n2+j] /= a;
+	  }
+	}
+
+	for (i = 0; i < n; i++)
+	  for (j = 0; j < n; j++)
+		dst->x[i*n+j] = r->x[i*n2+n+j];
+
+	return dst;
+}
+
+// Inspired from: http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis#Estimating_Regression_Models_Using_Least_Squares
+void multiple_reg_coeff(double *mx, double *my, int n, int k, double *coeff)
+{
+  	matrix X, Y;
+	X = mat_new(n,k);
+	X->x = mx;
+	Y = mat_new(n,1);
+	Y->x = my;
+
+    matrix mcoeff;
+	mcoeff = mat_mul(
+			mat_mul(
+				mat_inv(
+					mat_mul(transpose(X), X)),
+				transpose(X)),
+			Y);
+
+	for(int i=0; i<k; i++)
+		coeff[i] = mcoeff->x[i];
+
+	//mat_free(X);
+	//mat_free(Y);
+	//mat_free(mcoeff);
+}
+
+int test_multiple_regression()
+{
+	double da[] = {	1, 1,  1,   1,
+			2, 4,  8,  16,
+			3, 9, 27,  81,
+			4,16, 64, 256	};
+	double db[] = {     4.0,   -3.0,  4.0/3,
+			-13.0/3, 19.0/4, -7.0/3,
+			  3.0/2,   -2.0,  7.0/6,
+			 -1.0/6,  1.0/4, -1.0/6};
+
+	matrix_t a = { 4, 4, da }, b = { 4, 3, db };
+	matrix at;
+	at = transpose(&a);
+	matrix c = mat_mul(at, &b);
+
+	mat_show(&a), mat_show(at), mat_show(&b);
+	mat_show(c);
+	/* free(c) */
+	printf("\nInverse matrix:\n");
+
+
+	double dA[] = {	1, 2,  0,
+			-1, 1,  1,
+			1, 2, 3	};
+	matrix_t A = { 3, 3, dA };
+	mat_show(&A);
+	matrix Ainv;
+	Ainv = mat_inv(&A);
+	mat_show(Ainv);
+
+	// Multiple regression test: http://www.biddle.com/documents/bcg_comp_chapter4.pdf
+
+	double dX[] = {	1, 12, 32,
+			1, 14, 35,
+			1, 15, 45,
+			1, 16, 45,
+			1, 18, 50 };
+
+	double dY[] = {	350000, 399765, 429000, 435000, 433000};
+	int n = 5;
+	int k = 3;
+	matrix_t X= {5,k, dX};
+	matrix_t Y= {5,1, dY};
+	printf("\nMultiple regression:\n");
+	mat_show(&X);
+	mat_show(&Y);
+
+	matrix coeff;
+	coeff = mat_mul(
+			mat_mul(
+				mat_inv(
+					mat_mul(transpose(&X), &X)
+				        ),
+				transpose(&X)),
+			&Y);
+	mat_show(coeff);
+
+	double *results=NULL;
+	multiple_reg_coeff(dX, dY, n, k, results);
+	printf("\nFinal coefficients:\n");
+	for(int i=0; i<k; i++)
+	  printf("\tcoeff[%d]=%lf\n", i, results[i]);
+	return 0;
+
+}
+
+int _starpu_multiple_regression(struct starpu_perfmodel_history_list *ptr, double *coeff, unsigned ncoeff, unsigned nparameters, unsigned **combinations)
+{
+	unsigned n = find_list_size(ptr);
+	STARPU_ASSERT(n);
+
+	double *mx = (double *) malloc(ncoeff*n*sizeof(double));
+	STARPU_ASSERT(mx);
+
+	double *my = (double *) malloc(n*sizeof(double));
+	STARPU_ASSERT(my);
+
+	dump_multiple_regression_list(mx, my, ncoeff, nparameters, combinations, ptr);
+
+	multiple_reg_coeff(mx, my, n, ncoeff, coeff);
+
+	return 0;
+}

src/core/perfmodel/multiple_regression.h

@@ -0,0 +1,29 @@
+/* StarPU --- Runtime system for heterogeneous multicore architectures.
+ *
+ * Copyright (C) 2009, 2010  Université de Bordeaux
+ * Copyright (C) 2010  CNRS
+ *
+ * StarPU is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation; either version 2.1 of the License, or (at
+ * your option) any later version.
+ *
+ * StarPU 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 Lesser General Public License in COPYING.LGPL for more details.
+ */
+
+#ifndef __MULTIPLE_REGRESSION_H__
+#define __MULTIPLE_REGRESSION_H__
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <core/perfmodel/perfmodel.h>
+#include <starpu.h>
+
+int _starpu_multiple_regression(struct starpu_perfmodel_history_list *ptr, double *coeff, unsigned ncoeff, unsigned nparameters, unsigned **combinations);
+
+#endif // __MULTIPLE_REGRESSION_H__

src/core/perfmodel/perfmodel_history.c

@@ -31,6 +31,7 @@
 #include <core/workers.h>
 #include <datawizard/datawizard.h>
 #include <core/perfmodel/regression.h>
+#include <core/perfmodel/multiple_regression.h>
 #include <common/config.h>
 #include <starpu_parameters.h>
 #include <common/uthash.h>
@@ -288,25 +289,32 @@ static void dump_reg_model(FILE *f, struct starpu_perfmodel *model, int comb, in
 	 * Multiple Regression Model
 	 */
 
-	if (reg_model->ncoeff==0)
-		reg_model->ncoeff = model->ncombinations + 1;
-
-	reg_model->coeff = (double *) malloc(reg_model->ncoeff*sizeof(double));
 	if (model->type == STARPU_MULTIPLE_REGRESSION_BASED)
+	{
+		if (reg_model->ncoeff==0)
+			reg_model->ncoeff = model->ncombinations + 1;
+
+		reg_model->coeff = (double *) malloc(reg_model->ncoeff*sizeof(double));
 		_starpu_multiple_regression(per_arch_model->list, reg_model->coeff, reg_model->ncoeff, model->nparameters, model->combinations);
 
-	fprintf(f, "# n\tintercept");
-	for (int i=1; i < reg_model->ncoeff; i++){
-		fprintf(f, "\tc%d", i);
-	}
-	fprintf(f, "\n%u", reg_model->ncoeff);
-	for (int i=0; i < reg_model->ncoeff; i++){
-		fprintf(f, "\t%-15e", reg_model->coeff[i]);
+		fprintf(f, "# n\tintercept\t");
+		for (int i=0; i < model->ncombinations; i++)
+		{
+			if (model->parameters_names != NULL)
+				fprintf(f, "%s\t\t", model->parameters_names[i]);
+			else
+				fprintf(f, "c%d\t\t", i+1);
+		}
+
+		fprintf(f, "\n%u", reg_model->ncoeff);
+		for (int i=0; i < reg_model->ncoeff; i++)
+			fprintf(f, "\t%-15e", reg_model->coeff[i]);
+
 	}
 }
 #endif
 
-static void scan_reg_model(FILE *f, struct starpu_perfmodel_regression_model *reg_model)
+static void scan_reg_model(FILE *f, struct starpu_perfmodel_regression_model *reg_model, enum starpu_perfmodel_type model_type)
 {
 	int res;
 
@@ -351,22 +359,25 @@ static void scan_reg_model(FILE *f, struct starpu_perfmodel_regression_model *re
 	/*
 	 * Multiple Regression Model
 	 */
-	_starpu_drop_comments(f);
+	if (model_type == STARPU_MULTIPLE_REGRESSION_BASED)
+	{
+		_starpu_drop_comments(f);
 
-	// Read how many coefficients is there
-	res = fscanf(f, "%u", &reg_model->ncoeff);
-	STARPU_ASSERT_MSG(res == 1, "Incorrect performance model file");
+		// Read how many coefficients is there
+		res = fscanf(f, "%u", &reg_model->ncoeff);
+		STARPU_ASSERT_MSG(res == 1, "Incorrect performance model file");
 
-	reg_model->coeff = malloc(reg_model->ncoeff*sizeof(double));
+		reg_model->coeff = malloc(reg_model->ncoeff*sizeof(double));
 
-	unsigned multi_invalid = 0;
+		unsigned multi_invalid = 0;
 
-	for (int i=0; i < reg_model->ncoeff; i++){
-		res = _starpu_read_double(f, "%le", &reg_model->coeff[i]);
-		STARPU_ASSERT_MSG(res == 1, "Incorrect performance model file");
-		multi_invalid = (multi_invalid||isnan(reg_model->coeff[i]));
+		for (int i=0; i < reg_model->ncoeff; i++){
+			res = _starpu_read_double(f, "%le", &reg_model->coeff[i]);
+			STARPU_ASSERT_MSG(res == 1, "Incorrect performance model file");
+			multi_invalid = (multi_invalid||isnan(reg_model->coeff[i]));
+		}
+		reg_model->multi_valid = !multi_invalid;
 	}
-	reg_model->multi_valid = !multi_invalid;
 }
 
 
@@ -424,7 +435,7 @@ static void scan_history_entry(FILE *f, struct starpu_perfmodel_history_entry *e
 	}
 }
 
-static void parse_per_arch_model_file(FILE *f, struct starpu_perfmodel_per_arch *per_arch_model, unsigned scan_history)
+static void parse_per_arch_model_file(FILE *f, struct starpu_perfmodel_per_arch *per_arch_model, unsigned scan_history, enum starpu_perfmodel_type model_type)
 {
 	unsigned nentries;
 
@@ -433,7 +444,7 @@ static void parse_per_arch_model_file(FILE *f, struct starpu_perfmodel_per_arch
 	int res = fscanf(f, "%u\n", &nentries);
 	STARPU_ASSERT_MSG(res == 1, "Incorrect performance model file");
 
-	scan_reg_model(f, &per_arch_model->regression);
+	scan_reg_model(f, &per_arch_model->regression, model_type);
 
 	/* parse entries */
 	unsigned i;
@@ -492,7 +503,7 @@ static void parse_arch(FILE *f, struct starpu_perfmodel *model, unsigned scan_hi
 		{
 			struct starpu_perfmodel_per_arch *per_arch_model = &model->state->per_arch[comb][impl];
 			model->state->per_arch_is_set[comb][impl] = 1;
-			parse_per_arch_model_file(f, per_arch_model, scan_history);
+			parse_per_arch_model_file(f, per_arch_model, scan_history, model->type);
 		}
 	}
 	else
@@ -503,7 +514,7 @@ static void parse_arch(FILE *f, struct starpu_perfmodel *model, unsigned scan_hi
 	/* if the number of implementation is greater than STARPU_MAXIMPLEMENTATIONS
 	 * we skip the last implementation */
 	for (i = impl; i < nimpls; i++)
-		parse_per_arch_model_file(f, &dummy, 0);
+		parse_per_arch_model_file(f, &dummy, 0, model->type);
 
 }
 
@@ -1291,11 +1302,14 @@ double _starpu_multiple_regression_based_job_expected_perf(struct starpu_perfmod
 	reg_model = &model->state->per_arch[comb][nimpl].regression;
 
 	double parameter_value;
+	double *parameters;
+	parameters = (double *) malloc(model->nparameters*sizeof(double));
+	model->parameters(j->task, parameters);
 	expected_duration=reg_model->coeff[0];
 	for (int i=0; i < model->ncombinations; i++){
 		parameter_value=1.;
 		for (int j=0; j < model->nparameters; j++)
-			parameter_value *= pow(model->parameters[j],model->combinations[i][j]);
+			parameter_value *= pow(parameters[j],model->combinations[i][j]);
 
 		expected_duration += reg_model->coeff[i+1]*parameter_value;
 	}
@@ -1312,7 +1326,13 @@ docal:
 		model->benchmarking = 1;
 	}
 
-	return expected_duration*1000;//Time in milliseconds
+	// In the unlikely event that predicted duration is negative
+	// in case multiple linear regression is not so accurate
+	if (expected_duration < 0 )
+		expected_duration = 0.00001;
+
+	//Make sure that the injected time is in milliseconds
+	return expected_duration;
 }
 
 double _starpu_history_based_job_expected_perf(struct starpu_perfmodel *model, struct starpu_perfmodel_arch* arch, struct _starpu_job *j,unsigned nimpl)
@@ -1561,8 +1581,7 @@ void _starpu_update_perfmodel_history(struct _starpu_job *j, struct starpu_perfm
 			STARPU_ASSERT(entry);
 
 			entry->parameters = (double *) malloc(model->nparameters*sizeof(double));
-			for(int i=0; i < model->nparameters; i++)
-				entry->parameters[i] = model->parameters[i];
+			model->parameters(j->task, entry->parameters);
 			entry->duration = measured;
 
 			struct starpu_perfmodel_history_list *link;

src/core/perfmodel/regression.h

@@ -25,6 +25,5 @@
 #include <starpu.h>
 
 int _starpu_regression_non_linear_power(struct starpu_perfmodel_history_list *ptr, double *a, double *b, double *c);
-int _starpu_multiple_regression(struct starpu_perfmodel_history_list *ptr, double *coeff, unsigned ncoeff, unsigned nparameters, unsigned **combinations);
 
 #endif // __REGRESSION_H__ 

src/core/simgrid.c

@@ -326,14 +326,14 @@ void _starpu_simgrid_submit_job(int workerid, struct _starpu_job *j, struct star
 		return;
 
 	//LUKA ScalFMM injection
-	if (strcmp(_starpu_job_get_model_name(j),"p2p_cl_inout") == 0)
+	/*if (strcmp(_starpu_job_get_model_name(j),"p2p_cl_inout") == 0)
 	{
 	   length = 0.01582714*1000; //Time is in milliseconds
-	}
-	if (strcmp(_starpu_job_get_model_name(j),"p2p_cl_in") == 0)
+	}*/
+	/*if (strcmp(_starpu_job_get_model_name(j),"p2p_cl_in") == 0)
 	{
 	   length = 0.2636544*1000; //Time is in milliseconds
-	}
+	}*/
 	if (strcmp(_starpu_job_get_model_name(j),"p2m_cl") == 0)
 	{
 	   length = 1.474784*1000; //Time is in milliseconds