Welcome to MilkyWay@home

Source code improvement discussion


Advanced search

Message boards : Number crunching : Source code improvement discussion
Message board moderation

To post messages, you must log in.

1 · 2 · Next

AuthorMessage
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5392 - Posted: 9 Oct 2008, 5:04:50 UTC

Apologies on the length of this, but I felt it best to post some of this source code that is being said to be sloppy, since not many people know what is being talked about exactly.

The first example will be from the offending file that has the qgaus function, numericalIntegration.c:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define EPS 3.0e-11
#define PI 3.1415926535897932384626433832795028841971693993751

/*Gauss-Legendre quadrature taken from Numerical Recipes in C*/
void gaussLegendre(double x1, double x2, double x[], double w[], int n)
{
	int m, j, i;
	double z1,z,xm,xl,pp,p3,p2,p1;

	m = (n+1)/2;
	xm = 0.5*(x2+x1);
	xl = 0.5*(x2-x1);
	
	//fprintf(stderr, "m = %d: xm = %g: xl = %gn", m, xm, xl);
	for (i=1; i<=m; i++) {
		//fprintf(stderr, "starting iteration %d of outer loopn", i);
		z = cos(PI*(i-0.25)/(n+0.5));
		do {
			p1 = 1.0;
			p2 = 0.0;
			for (j = 1; j <= n; j++) {
				//fprintf(stderr, "starting iteration %d of inner loopn", i);
				p3 = p2;
				p2 = p1;
				p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp = n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z = z1 - p1/pp;
			//fprintf(stderr, "z-z1 = %gn", fabs(z-z1));
		} while (fabs(z-z1) > EPS);
		
		x[i-1] = xm - xl *z;
		x[n-i] = xm + xl *z;
		w[i-1] = 2.0*xl/((1.0-z*z)*pp*pp);
		w[n-i] = w[i-1];
	}
}

double qgaus(double (*func)(double), double a, double b, int n) {
	int j;
	double xr, xm, dx, s;
	double *x;
	double *w;

	x = (double*)malloc(sizeof(double)*n);
	w = (double*)malloc(sizeof(double)*n);

	gaussLegendre(-1.0, 1.0, x, w, n);

	xm = 0.5*(b+a);
	xr = 0.5*(b-a);
	s = 0;
	for (j = 0; j < n; j++) {
		dx = xr*x[j];
		s += w[j]*((*func)(xm+dx));
	}

	free(x);
	free(w);
	return s *= xr;
}


Immediately what I see is not the math, but the fact that the variables have no immediate meaning to someone not intimately involved in the process. Like I said, I have been involved in business logic and not a lot of math for quite some time. I also am not real familiar with C. I did a little C++ and knowing Java helps. The other major thing missing are comments.

Next is the code from integral_function.c. When you look at it, notice who it says wrote it and notice the complete difference in style...

/*
 *  integral_function.c
 *  Astronomy
 *
 *  Created by Travis Desell on 2/21/07.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

/****
        *       BOINC includes
*****/

#ifdef _WIN32
#include "boinc_win.h"
#else
#include "config.h"
#endif

#ifndef _WIN32
#include <cstdio>
#include <cctype>
#include <cstring>
#include <cstdlib>
#include <csignal>
#endif

//#define BOINC_APP_GRAPHICS

#ifdef BOINC_APP_GRAPHICS
#include "graphics_api.h"
#include "graphics_lib.h"
#endif

#include "diagnostics.h"
#include "str_util.h"
#include "util.h"
#include "filesys.h"
#include "boinc_api.h"
#include "mfile.h"

using std::string;

#define INTEGRAL_CHECKPOINT_FILE "astronomy_integral_checkpoint"

/****
	*	Astronomy includes
*****/
#include "integral_function.h"
#include "parameters.h"
#include "volume.h"
#include "probability.h"
#include "stCoords.h"
#include "atSurveyGeometry.h"
#include "star_points.h"

#ifndef _WIN32
#define pi M_PI
#else
#define pi 3.14159265358979323846
#endif
#define deg (180.0/pi)

bool convolve = false;
double integral_time;

double background_integral;
double *stream_integrals;

int mu_step_current, nu_step_current, r_step_current;

double get_background_integral() {
	return background_integral;
}

double* get_stream_integrals() {
	return stream_integrals;
}

double get_integral_time() {
	return integral_time;
}

void load_integral_checkpoint() {
	FILE* checkpoint;
	char chkpt_path[512];

	fprintf(stderr,"APP: astronomy reading integral checkpoint filen");

	// See if there's a valid checkpoint file.
	// If so seek input file
	boinc_resolve_filename(INTEGRAL_CHECKPOINT_FILE, chkpt_path, sizeof(chkpt_path));
	checkpoint = boinc_fopen(chkpt_path, "r");
	if (checkpoint) {
		//double likelihood;
		fscanf(checkpoint, "mu_step_current: %d, nu_step_current: %d, r_step_current: %d, background_integral: %lfn", &mu_step_current, &nu_step_current, &r_step_current, &background_integral);
		stream_integrals = (double*)malloc(sizeof(double) * get_number_streams());
		read_double_array(checkpoint, "stream_integrals: ", get_number_streams(), &stream_integrals);
		fclose(checkpoint);
	} else {
		mu_step_current = 0;
		nu_step_current = 0;
		r_step_current = 0;
		background_integral = 0;
		stream_integrals = (double*)malloc(sizeof(double) * get_number_streams());
		for (int i = 0; i < get_number_streams(); i++) {
			stream_integrals[i] = 0;
		}
	}

	fprintf(stderr,"APP: astronomy read integral checkpoint finishedn");
}

int do_integral_checkpoint() {
	string resolved_name;
	int retval;

	FILE* checkpoint = fopen("temp", "w");

	if (checkpoint) {
		fprintf(checkpoint, "mu_step_current: %d, nu_step_current: %d, r_step_current: %d, background_integral: %lfn", mu_step_current, nu_step_current, r_step_current, background_integral); 
		print_double_array(checkpoint, "stream_integrals: ", get_number_streams(), stream_integrals);
		fclose(checkpoint);
		//return 0;
	} else {
		fprintf(stderr, "APP: astronomy ERROR GENERATING CHECKPOINT FILEn");
		return 1;
	}

	fprintf(stderr,"APP: astronomy integral checkpointingn");

	boinc_resolve_filename_s(INTEGRAL_CHECKPOINT_FILE, resolved_name);
	retval = boinc_rename("temp",resolved_name.c_str());
	if(retval) return retval;

	fprintf(stderr,"APP: astronomy integral checkpoint donen");
	return 0;
}

void calculate_integrals() {
	int wedge;
	double background_weight;
	double *background_parameters;
	double *stream_weights;
	double **stream_parameters;

	double nu, nu_min, nu_max, nu_step, nu_steps;
	double mu, mu_min, mu_max, mu_step, mu_steps;
	double r, r_min, r_max, r_step, r_steps, log_r, next_r;

	int num_calls;
	double total_volume, V;
	double ir, id;
	double ra, dec, point0, point1;
	double integral_point[3], xyz[3];

	int retval;
	double f;
	bool first_run = true;

	wedge = get_wedge();
	background_weight = get_background_weight();
	background_parameters = get_background_parameters();
	stream_weights = get_stream_weights();
	stream_parameters = get_stream_parameters();

	nu_min = get_nu_min();
	nu_max = get_nu_max();
	nu_step = get_nu_step();
	nu_steps = get_nu_steps();
	
	mu_min = get_mu_min();
	mu_max = get_mu_max();
	mu_step = get_mu_step();
	mu_steps = get_mu_steps();
	
	r_min = get_r_min();
	r_max = get_r_max();
	r_step = get_r_step();
	r_steps = get_r_steps();

	background_integral = 0;
	num_calls = 0;
	V = 0;
	total_volume = 0;

	//printf("r_min: %lf, r_max: %lf, r_step: %lf, r_steps: %lfn", r_min, r_max, r_step, r_steps);
	//printf("mu_min: %lf, mu_max: %lf, mu_step: %lf, mu_steps: %lfn", mu_min, mu_max, mu_step, mu_steps);
	//printf("nu_min: %lf, nu_max: %lf, nu_step: %lf, nu_steps: %lfn", nu_min, nu_max, nu_step, nu_steps);
	//printf("wedge: %dn", wedge);

	//printf("background_weight: %lfn", background_weight);
	//print_double_array(stdout, "background_parameters: ", 4, background_parameters);
	//printf("stream_weight: %lfn", stream_weights[0]);
	//print_double_array(stdout, "stream_parameters: ", 5, stream_parameters[0]);
	//printf("n");

	load_integral_checkpoint();

	for ( ; mu_step_current < mu_steps; mu_step_current++) {
		mu = mu_min + mu_step_current*mu_step;

		if(!first_run) nu_step_current=0;
		for ( ; nu_step_current < nu_steps; nu_step_current++) {
			nu = nu_min + nu_step_current*nu_step;

			if(!first_run) r_step_current=0;
			for ( ; r_step_current < r_steps; r_step_current++) {
				log_r = r_min + r_step_current*r_step;

				r = pow(10.0, (log_r-14.2)/5.0);
				next_r = pow(10.0, (log_r+r_step-14.2)/5.0);

				if (wedge > 0) {
					ir = (pow(next_r,3.0) - pow(r, 3.0))/3.0;
					id = cos((90 - nu - nu_step)/deg) - cos((90 - nu)/deg);

					V = ir * id * mu_step / deg;

					ra = 0.0;
					dec = 0.0;
					point0 = 0.0;
					point1 = 0.0;

					atGCToEq(mu + 0.5 * mu_step, nu + 0.5 * nu_step, &ra, &dec, get_node(), wedge_incl(wedge));
					atEqToGal(ra, dec, &point0, &point1);
					integral_point[0] = point0;
					integral_point[1] = point1;
					integral_point[2] = (next_r+r)/2.0;

					total_volume += V;
				} else {
					V = mu_step * nu_step * r_step;
					xyz[0] = mu + 0.5*mu_step;
					xyz[1] = nu + 0.5*nu_step;
					xyz[2] = r + 0.5*r_step;
					xyz2lbr(xyz, integral_point);
				}

				double bg_prob = 0.0;
				double st_prob = 0.0;

				if (convolve) {
//					bg_prob = Probability.stPbxConvolved(integral_point, background_parameters);
//					st_prob = Probability.stPsgConvolved(integral_point, stream_parameters[0], wedge);
				} else {
					bg_prob = stPbx(integral_point, background_parameters);
					st_prob = stPsg(integral_point, stream_parameters[0]);
				}

				background_integral += bg_prob * V;
				stream_integrals[0] += st_prob * V;
				first_run = false;
				num_calls++;

				/*if (num_calls%100000 == 0) {
					fprintf(stderr, "integral point: %lf %lf %lfn", integral_point[0], integral_point[1], integral_point[2]);
					fprintf(stderr, "ir: %lf, id: %lf, deg: %lfn", ir, id, deg);
					fprintf(stderr, "V: %lf, bg_prob: %lf, st_prob: %lfn", V, bg_prob, st_prob);
					fprintf(stderr, "bg_prob: %lf, st_prob: %lf, total_volume: %lfn", bg_prob, st_prob, total_volume);
					fprintf(stderr, "completed %d calls. background_integral: %lf, stream_integral[0]: %lfn", num_calls, background_integral, stream_integral[0]);
					}*/

				if(boinc_time_to_checkpoint()) {
				        retval = do_integral_checkpoint();

					if (retval) {
					          fprintf(stderr,"APP: astronomy checkpoint failed %dn",retval);
						  exit(retval);
					}
					boinc_checkpoint_completed();
				}

				f = mu_step_current*nu_steps*r_steps + nu_step_current*r_steps + r_step_current;
				f /= mu_steps*nu_steps*r_steps + get_number_streams()*(get_number_star_points()+.001)/100;
				//if(num_calls%100000 == 0) printf("%f percent donen",f*100);
				boinc_fraction_done(f);
			}
			r_step_current = 0;
		}
		nu_step_current = 0;
	}

	do_integral_checkpoint();
}


Meaningful variable names, like ra and dec, meaning Right Ascension and Declination. Even if you don't comprehend the math, you at least have some general idea of the basic tasks that are being performed. I don't know if they are still working on that particular section of code or not. What is missing from that example is a changelog...

Unless the first example was really code that Travis wrote and just doesn't want to own up to, which means I would be calling him a liar if I said that, then I'm inclined to believe what Travis posted, that they didn't write it and are working on other parts of the application.

Brian
ID: 5392 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5395 - Posted: 9 Oct 2008, 9:38:15 UTC - in response to Message 5392.  
Last modified: 9 Oct 2008, 9:51:34 UTC

The first example will be from the offending file that has the qgaus function, numericalIntegration.c:

[..]

Immediately what I see is not the math, but the fact that the variables have no immediate meaning to someone not intimately involved in the process. Like I said, I have been involved in business logic and not a lot of math for quite some time. I also am not real familiar with C. I did a little C++ and knowing Java helps. The other major thing missing are comments.

Next is the code from integral_function.c. When you look at it, notice who it says wrote it and notice the complete difference in style...

[..]

Meaningful variable names, like ra and dec, meaning Right Ascension and Declination. Even if you don't comprehend the math, you at least have some general idea of the basic tasks that are being performed. I don't know if they are still working on that particular section of code or not. What is missing from that example is a changelog...

Unless the first example was really code that Travis wrote and just doesn't want to own up to, which means I would be calling him a liar if I said that, then I'm inclined to believe what Travis posted, that they didn't write it and are working on other parts of the application.

Brian

Come on Brian, you could do better.

The names in qgaus are not meaningful? It is just doing a one dimensional integration of a function (google for Gauss-Legendre).
How one calls the coordinate in one dimension? Yes, simply x. If you deal with an integration (or differentiation) and have the slightests mathematical background, you should know that the differential x is commonly abbreviated as dx. xm = (a+b)* 0.5 is simply the mean value of a and b and the actual x value we are looking for (x-mean, or short xm). If you now finished googling for Gauss-Legendre quadrature, you should know, that w stands for weight. s is of course the sum, and finally n is just the number of points we are using to approximate the integral. Was that so difficult? I think not more than to figure out what ra, dec, GC, Eq, lbr, mu, or nu actually mean.

I think, the naming of the variables isn't the big problem there. But one should actually look like, what values these "variables" could have. One finds quite fast, that most of them are actually constant (that means the same for the whole runtime of the app). Especially is n (numpoints) constant, that means you can declare x and w on filescope and initialize once and you are done. You can spare the call to GaussLegendre every time, which speeds up the whole thing tremendously.
Further things to note are that b-a and therefore dx is constant, xm = (a+b)*0.5 is known before calling qgaus and a and b are only calculated immediately before calling qgaus. For what? Why not passing just xm?

By the way, the version of "integral_function.c" you showed is not functional (the convolution stuff is commented out, that part is actually doing the work, as convolve is 30 in current WUs). If it is a part of the new version of the app, I have to say, that not much have improved. In the current source it is named "evaluation.c" and almost identical.

Edit:
I've just seen in the other thread where you got this code from. It's an outdated version. The current source is only 7 and a half month old, not a full year.
ID: 5395 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5397 - Posted: 9 Oct 2008, 11:40:56 UTC - in response to Message 5395.  
Last modified: 9 Oct 2008, 11:48:05 UTC

The first example will be from the offending file that has the qgaus function, numericalIntegration.c:

[..]

Immediately what I see is not the math, but the fact that the variables have no immediate meaning to someone not intimately involved in the process. Like I said, I have been involved in business logic and not a lot of math for quite some time. I also am not real familiar with C. I did a little C++ and knowing Java helps. The other major thing missing are comments.

Next is the code from integral_function.c. When you look at it, notice who it says wrote it and notice the complete difference in style...

[..]

Meaningful variable names, like ra and dec, meaning Right Ascension and Declination. Even if you don't comprehend the math, you at least have some general idea of the basic tasks that are being performed. I don't know if they are still working on that particular section of code or not. What is missing from that example is a changelog...

Unless the first example was really code that Travis wrote and just doesn't want to own up to, which means I would be calling him a liar if I said that, then I'm inclined to believe what Travis posted, that they didn't write it and are working on other parts of the application.

Brian

Come on Brian, you could do better.


Perhaps, but as I said, I haven't been involved in anything more than simple ("business") math in a long time. I think the last Calculus class I took was in 1993, but I was barely paying attention in that class, so the last time I think I was actually paying attention to higher level math was perhaps 1991. I don't know how old you are, but I'm 37. Anyone will lose the mental skills that one had when they were younger if those skills are not in use.

At any rate, the point I was trying to make is you can see two totally different coding styles in the two examples. Travis has stated that he and Nate didn't write it and are not going through that code right now. I, personally, believe them.

What you can do to help me and others believe you is for you to perhaps take this code and explain, in detail, how it is inefficient and what would be efficient. No "you could do better" or make mention about "anyone with the slightest mathematical background should know". See, there is an underlying level of mistrust that what you have done is completely honest. Not by me, but by those people who are raising issues of things not being fair. You'd improve your position with everyone if you took some time to explain things in detail, lest you start seeming like Francois ("Who?") over at SETI... Those of us who were around for all the marketing that he did will know what I'm talking about...

I gotta go...

Brian

Edit: I see you did some explaining, but I'm honestly in the middle of getting ready to leave, so I haven't had time to look at what you said and think about it. I'll do that sometime today...
ID: 5397 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Opteron

Send message
Joined: 1 Oct 08
Posts: 5
Credit: 1,892
RAC: 0
1 credit badge10 year member badge
Message 5398 - Posted: 9 Oct 2008, 12:08:41 UTC - in response to Message 5397.  

What you can do to help me and others believe you is for you to perhaps take this code and explain, in detail, how it is inefficient and what would be efficient.

Well he has already one example (probably the worst) in his profile:
To cite the Milkyway source code:

sint = (zp * sinl) / d;
xyz[1] = d * sint;


written in a more easy way with shorter variable names this would be:

x = a * b / d;
return_value = x * d;

It should be obvious to most people that the division and multiplication with the same parameter "d" is .. well ... totally unnecessary.

The shortest way to write the above statement in a better way would be:
return_value = a * b;

I believe him if he says that there is more like this, because this example tells already a lot about the general code quality.

Have a nice day

Opteron
ID: 5398 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5399 - Posted: 9 Oct 2008, 12:23:10 UTC - in response to Message 5397.  

Perhaps, but as I said, I haven't been involved in anything more than simple ("business") math in a long time. I think the last Calculus class I took was in 1993, but I was barely paying attention in that class, so the last time I think I was actually paying attention to higher level math was perhaps 1991. I don't know how old you are, but I'm 37. Anyone will lose the mental skills that one had when they were younger if those skills are not in use.

At any rate, the point I was trying to make is you can see two totally different coding styles in the two examples. Travis has stated that he and Nate didn't write it and are not going through that code right now. I, personally, believe them.

What you can do to help me and others believe you is for you to perhaps take this code and explain, in detail, how it is inefficient and what would be efficient. No "you could do better" or make mention about "anyone with the slightest mathematical background should know". See, there is an underlying level of mistrust that what you have done is completely honest. Not by me, but by those people who are raising issues of things not being fair. You'd improve your position with everyone if you took some time to explain things in detail, lest you start seeming like Francois ("Who?") over at SETI... Those of us who were around for all the marketing that he did will know what I'm talking about...

I gotta go...

Brian

Edit: I see you did some explaining, but I'm honestly in the middle of getting ready to leave, so I haven't had time to look at what you said and think about it. I'll do that sometime today...

If you want to know my age, I'm 30. My partner at this account (who has also done the "optimizations" independently from me) is slightly older. And of course one forgets things one haven't done for quite some time. That was also the reason I explicitly stated, that I'm not a good coder and haven't done any programming for some years. Still it was easy to spot and correct some very inefficient things in a day or so (after downloading the source code, after my normal work on one evening). And remember, I am not familiar with the stuff done in the code and have of course not written it.

The comment "you could do better" was referring to your comparison of the two c files. Frankly, I don't see too much difference in style there. And when you could figure out, what ra, dec, GC, Eq, lbr, mu, or nu actually mean, why you were not able to do it in the qgaus function? ra for right ascension or dec for declination is understandable on the same level as w for weight or dx for differential x. No difference here.

Concerning my explanations, I already said that I don't want to take Travis, Nate and Dave to school. My explanations should really be adequate and more than sufficient for a computer scientist. Remember, we don't want to do here a beginners course in programming. That would be somehow an insult to the project staff. I can't believe they don't have the skills. They just have to take a day or two and do the changes. We are waiting for half a year already.

And frankly, I don't think I sound like Francois (who?). I have already discussed with him over at aceshardware. He has quite some knowledge, but he also takes sometimes extreme positions, which are just not logical and fail to explain it even to some guys with a comparable knowledge. Here, it is different. Every somehow capable programmer I know of has the same opinion about the MW source.
ID: 5399 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5400 - Posted: 9 Oct 2008, 12:30:45 UTC - in response to Message 5398.  

Well he has already one example (probably the worst) in his profile:
To cite the Milkyway source code:

sint = (zp * sinl) / d;
xyz[1] = d * sint;


written in a more easy way with shorter variable names this would be:

x = a * b / d;
return_value = x * d;

It should be obvious to most people that the division and multiplication with the same parameter "d" is .. well ... totally unnecessary.

The shortest way to write the above statement in a better way would be:
return_value = a * b;

I believe him if he says that there is more like this, because this example tells already a lot about the general code quality.

Have a nice day

Opteron

Well, it is not the worst example from the performance standpoint, it is only one of the absurdest. The qgaus function cost far more performance than this division.
The main problem is really that a lot of the stuff is calculated not only twice, but millions and millions times. So more like the example with the nested for-loops. Everything you can shift from the innermost loop to an outer loop will help a lot. And this is a very easy thing to understand and shouldn't need much explanation to a computer scientist.
ID: 5400 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5401 - Posted: 9 Oct 2008, 13:47:28 UTC - in response to Message 5399.  
Last modified: 9 Oct 2008, 13:50:57 UTC

That was also the reason I explicitly stated, that I'm not a good coder and haven't done any programming for some years. Still it was easy to spot and correct some very inefficient things in a day or so (after downloading the source code, after my normal work on one evening). And remember, I am not familiar with the stuff done in the code and have of course not written it.


...and perhaps I might spot stuff easily too, but I haven't been looking at it. That is what they are saying to you: they are not looking at it right now and have other priorities that have been deemed by either themselves or by a project team / leader as higher priorities than this...

You mention that the code I posted is older. If you'd like, feel free to PM me and I'll give you my email address where you can send the unmodified base code that you started from to me and I'll take some time this weekend and look at it. Can't promise I'll understand or catch what you see as inefficient, but I'll take a look. It might change my opinion on the matter...


The comment "you could do better" was referring to your comparison of the two c files. Frankly, I don't see too much difference in style there.


You apparently don't have to do much maintaining of code. At the job before the one I have now, we had to count lines of code for some audit. We didn't manually do it, so it was an estimate, but I want to say it was close to 1 million lines of code for just one portion of the total software package. Readability is extremely nice.

And when you could figure out, what ra, dec, GC, Eq, lbr, mu, or nu actually mean, why you were not able to do it in the qgaus function?


...because as a child, I had a reflector and a refractor telescope. As a teen, I went up into the "mountains" (in quotes because they're small hills in compared to other mountains) and was the first person out of a group of around 200 people to spot Halley's Comet... That said, however, GC, Eq, lbr, mu, nu, and other Physics-related variables, I'm extremely rusty on, same thing for the stuff like differential or momenta.

To illustrate why I still believe the code in numericalIntegration is not documented properly first and foremost, I offer the following code that was obtained from here
Edit: Looking at that, I also spot some silly use of block commenting, but that's at least not harmful, since they are comments...but the block comment doesn't need to be ended on individual lines. If they wanted to do that, then just doing // on each line would've accomplished the same thing...

/*******************************************************************************
Gauss-Legendre integration function, gauleg, from "Numerical Recipes in C"
(Cambridge Univ. Press) by W.H. Press, S.A. Teukolsky, W.T. Vetterling, and
B.P. Flannery
*******************************************************************************/

#include <stdlib.h>
#include <math.h>

#define NR_END 1
#define EPS 3.0e-11 /* EPS is the relative precision. */

double *dvector(int nl, int nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;
	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
	return v-nl+NR_END;
}

/******************************************************************************/
void gauleg(double x1, double x2, double x[], double w[], int n)
/*******************************************************************************
Given the lower and upper limits of integration x1 and x2, and given n, this
routine returns arrays x[1..n] and w[1..n] of length n, containing the abscissas
and weights of the Gauss-Legendre n-point quadrature formula.
*******************************************************************************/
{
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
	m=(n+1)/2; /* The roots are symmetric, so we only find half of them. */
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
	for (i=1;i<=m;i++) { /* Loop over the desired roots. */
		z=cos(3.141592654*(i-0.25)/(n+0.5));
		/* Starting with the above approximation to the ith root, we enter */
		/* the main loop of refinement by Newton's method.                 */
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) { /* Recurrence to get Legendre polynomial. */
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			/* p1 is now the desired Legendre polynomial. We next compute */
			/* pp, its derivative, by a standard relation involving also  */
			/* p2, the polynomial of one lower order.                     */
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp; /* Newton's method. */
		} while (fabs(z-z1) > EPS);
		x[i]=xm-xl*z;      /* Scale the root to the desired interval, */
		x[n+1-i]=xm+xl*z;  /* and put in its symmetric counterpart.   */
		w[i]=2.0*xl/((1.0-z*z)*pp*pp); /* Compute the weight             */
		w[n+1-i]=w[i];                 /* and its symmetric counterpart. */
	}
}


I have not compared the two to know which one might be more efficient or if they are both the same, just posting it to illustrate the higher level of documentation. Not sure about qgaus, but there could be "stock" code from a book or other online resource out there for that as well. Don't have time to look now though...

I can't believe they don't have the skills. They just have to take a day or two and do the changes. We are waiting for half a year already.


What is so difficult for you to understand when they say "we aren't working on that right now, if you'd like to help us out, please, help us out"?????? Has the thought occurred to you that maybe they would see the same things if they did look at it? They are stating that they have other priorities and have asked for and encouraged your assistance. You return that by essentially calling them "stupid" or "lazy", although somewhat politely...


Every somehow capable programmer I know of has the same opinion about the MW source.


Be that as it may, they have now repeatedly said that they are working on other things and have asked for your help. The choice is yours as to if you want to help or not.
ID: 5401 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5402 - Posted: 9 Oct 2008, 14:34:55 UTC - in response to Message 5401.  
Last modified: 9 Oct 2008, 14:37:15 UTC

To illustrate why I still believe the code in numericalIntegration is not documented properly first and foremost, I offer the following code that was obtained from here
Edit: Looking at that, I also spot some silly use of block commenting, but that's at least not harmful, since they are comments...but the block comment doesn't need to be ended on individual lines. If they wanted to do that, then just doing // on each line would've accomplished the same thing...

I agree that the code in numeicalIntegration.c is not documented properly. But you made the point that you think the other file (evaluation.c or integral_function.c) is better. I don't see an abundance of comments there either.
And concerning the block comments, as the code is written for pure ANSI C, the comments are fine. Originally, C did not accept commenting with double slashes ;)

I have not compared the two to know which one might be more efficient or if they are both the same, just posting it to illustrate the higher level of documentation.

It's the same aside from the comments. It is just copied from the book "Numerical Recipes in C".
What is so difficult for you to understand when they say "we aren't working on that right now, if you'd like to help us out, please, help us out"?????? Has the thought occurred to you that maybe they would see the same things if they did look at it? They are stating that they have other priorities and have asked for and encouraged your assistance. You return that by essentially calling them "stupid" or "lazy", although somewhat politely...

But on what are they working the last 6 months? Tuning the genetic algorithm producing the WUs? This part is not much more complex than the application (you can look at it too, I will send you a PM). Considering the overall importance of the app (over 99% of the total computation time is spent there, not in the genetic algorithm stuff), that may be just wrong priorities? Considering the low amount of work necessary to get massive speedups of the app (together with a cleaner and easier to maintain coding), it is simply not understandable what they wait for. Simply speaking, a waste of computing power of this dimension by this project can not be justified, not after months of standstill. We are not talking about a project just started, it runs already for quite some time.

PS:
The code tags suck, all the indentation is lost, hampering the readability. The pre tags are better for small code snippets, even if the text gets double spaced.
ID: 5402 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5405 - Posted: 9 Oct 2008, 15:10:45 UTC - in response to Message 5402.  


I have not compared the two to know which one might be more efficient or if they are both the same, just posting it to illustrate the higher level of documentation.

It's the same aside from the comments. It is just copied from the book "Numerical Recipes in C".


Am I to assume that you don't have any issues with that portion and only the qgaus code? Not asking as a "trick question". I haven't looked over it and, like I said, math is just not my strong suit anymore. The most I've had to do in at least 10 years is multiplication and division. When you're dealing with money (retail), you just aren't having to do curve fitting and other fancy stuff...

Considering the overall importance of the app (over 99% of the total computation time is spent there, not in the genetic algorithm stuff), that may be just wrong priorities?


I'm assuming you ran a profiler to get that? Since the system I use and would be comfortable mucking around on is an AMD, I can't use VTune, or at least I remember trying to use it and not having much luck...

As for pre vs. code, meh...now you're complaining about the formatting too? ;-P
ID: 5405 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5406 - Posted: 9 Oct 2008, 15:16:45 UTC - in response to Message 5402.  


I have not compared the two to know which one might be more efficient or if they are both the same, just posting it to illustrate the higher level of documentation.

It's the same aside from the comments. It is just copied from the book "Numerical Recipes in C".


Am I to assume that you don't have any issues with that portion and only the qgaus code? Not asking as a "trick question". I haven't looked over it and, like I said, math is just not my strong suit anymore. The most I've had to do in at least 10 years is multiplication and division. When you're dealing with money (retail), you just aren't having to do curve fitting and other fancy stuff...

Considering the overall importance of the app (over 99% of the total computation time is spent there, not in the genetic algorithm stuff), that may be just wrong priorities?


I'm assuming you ran a profiler to get that? Since the system I use and would be comfortable mucking around on is an AMD, I can't use VTune, or at least I remember trying to use it and not having much luck...

As for pre vs. code, meh...now you're complaining about the formatting too? ;-P
ID: 5406 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5407 - Posted: 9 Oct 2008, 15:35:36 UTC - in response to Message 5405.  
Last modified: 9 Oct 2008, 15:36:39 UTC

Am I to assume that you don't have any issues with that portion and only the qgaus code? Not asking as a "trick question". I haven't looked over it and, like I said, math is just not my strong suit anymore. The most I've had to do in at least 10 years is multiplication and division. When you're dealing with money (retail), you just aren't having to do curve fitting and other fancy stuff...

The GaussLegendre function is doing what it should do. The "Numerical Recipes" are actually quite a good starting point.

The problem is it gets called about 180 million times for a 260cr WU when in fact it only needs to be called exactly once.

I'm assuming you ran a profiler to get that? Since the system I use and would be comfortable mucking around on is an AMD, I can't use VTune, or at least I remember trying to use it and not having much luck...

No, just an estimate as it runs on the server and I know that a genetic algorithm is not that computationally intensive. Actually I think it is even more like 99.9+ percent.
I have not used VTune at all (I guess it would not help much at the state of the original code, but I have never seen it yet). As I said, my version is still unoptimized. I just cutted out things like unnecessary function calls and calculations. The only thing I used is the profiler integrated in VC98 (quite outdated, I know, but I have not done such stuff for years) to get an idea where the application spends all the time. And my box is an AthlonXP ;)

As for pre vs. code, meh...now you're complaining about the formatting too? ;-P

Consider that as a general remark ;) I think the indentation is really important for the readability.
ID: 5407 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5412 - Posted: 9 Oct 2008, 16:32:38 UTC - in response to Message 5407.  

Am I to assume that you don't have any issues with that portion and only the qgaus code? Not asking as a "trick question". I haven't looked over it and, like I said, math is just not my strong suit anymore. The most I've had to do in at least 10 years is multiplication and division. When you're dealing with money (retail), you just aren't having to do curve fitting and other fancy stuff...

The GaussLegendre function is doing what it should do. The "Numerical Recipes" are actually quite a good starting point.

The problem is it gets called about 180 million times for a 260cr WU when in fact it only needs to be called exactly once.

I'm assuming you ran a profiler to get that? Since the system I use and would be comfortable mucking around on is an AMD, I can't use VTune, or at least I remember trying to use it and not having much luck...

No, just an estimate as it runs on the server and I know that a genetic algorithm is not that computationally intensive. Actually I think it is even more like 99.9+ percent.
I have not used VTune at all (I guess it would not help much at the state of the original code, but I have never seen it yet). As I said, my version is still unoptimized. I just cutted out things like unnecessary function calls and calculations. The only thing I used is the profiler integrated in VC98 (quite outdated, I know, but I have not done such stuff for years) to get an idea where the application spends all the time. And my box is an AthlonXP ;)


At one point I had Visual Studio 6 Enterprise Edition on my old computer, but that system got mothballed when I got my AMD system and I have never installed it again. I need to check into getting Visual Studio 2005 or 2008. I do have the CDs for a beta version of 2005, but I think best not to go with a beta.

As for the code you posted in your profile:

/* Convert sun-centered lbr into galactic xyz coordinates. */
void lbr2xyz(const double* lbr, double* xyz) {
    double r0, sinb, sinl, cosb, cosl, sint, cost;
    double zp, d;
    
    r0 = 8.5;
    sinb = sin(lbr[1] / deg);
    sinl = sin(lbr[0] / deg);
    cosb = cos(lbr[1] / deg);
    cosl = cos(lbr[0] / deg);

    xyz[2] = lbr[2] * sinb;
    zp = lbr[2] * cosb;
    d = sqrt( r0 * r0 + zp * zp - 2 * r0 * zp *cosl );

    sint = (zp * sinl) / d;
    cost = (zp * zp - r0 * r0 - d * d) / (2 * d * r0);

    xyz[0] = d * cost;
    xyz[1] = d * sint;
}


The sine and cosine variables are not references / pointers, so there's no need in doing some things with them, and more efficient to do all math and assignment into the xyz[] array in one step vs. multiple steps...unless of course you do think that there will be a need for individual sine and cosine values down the road, but of course you could then create the variables and store that information for reuse in the same function.

I'm assuming that's similar to what you're finding elsewhere?
ID: 5412 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5415 - Posted: 9 Oct 2008, 17:40:23 UTC - in response to Message 5412.  

As for the code you posted in your profile:

[..]

The sine and cosine variables are not references / pointers, so there's no need in doing some things with them, and more efficient to do all math and assignment into the xyz[] array in one step vs. multiple steps...unless of course you do think that there will be a need for individual sine and cosine values down the road, but of course you could then create the variables and store that information for reuse in the same function.

I'm assuming that's similar to what you're finding elsewhere?

Oh, you have done the first part of the challenge for the project, you located the lines! Yes, that's the function.
I won't go into details, as this is a task I want to be done by the project staff. I will only say yes or no. But as a reminder, the fastest function is the function that does not get called. So the complete solution wouldn't only be a new version of this function (it seems you think in the right direction) but also a significantly reduced count of calls. In the current state it gets called more than 8 billion times for a 260cr WU. I want to see a significantly smaller number for that. That's actually the harder part of the challenge, I guess.

And yes, that is a typical example, also the amount of unnecessary calls.
ID: 5415 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5423 - Posted: 9 Oct 2008, 20:19:12 UTC - in response to Message 5415.  
Last modified: 9 Oct 2008, 20:29:36 UTC

So the complete solution wouldn't only be a new version of this function (it seems you think in the right direction) but also a significantly reduced count of calls. In the current state it gets called more than 8 billion times for a 260cr WU. I want to see a significantly smaller number for that. That's actually the harder part of the challenge, I guess.


The qgaus call isn't in a loop, but the function that it is a part of could be. Chasing up, I'm finding a set of loops with a struct involved. I'm guessing around without the benefit of being able to debug and step through...

Also, I do notice a difference between what I posted for gaussLegendre (gauleg) and what the project has after the while - the array positions. Again, without the capability to step it through, I don't know if it matters, but it appears to matter.

Edit: BTW, you are stating in bold two different figures, one of "180 million times", then another of "8 billion times". Without the benefit of being able to run the code in a debugger or have a profiler profile it, I can't definitively state that I think you're exaggerating things for dramatic effect, but clearly an approximate 44 times increase in the number of iterations that are deemed "excessive", particularly at that high of a number, seems unlikely. If it really goes through there "more than it should", then expressing that as either "more than it should" or with a number closer to what it really does would make your claim more believable... Over-dramatizing isn't helpful. Now if it really does loop 8 billion times, then yeah, there's a severe problem.
ID: 5423 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Milksop at try

Send message
Joined: 1 Oct 08
Posts: 106
Credit: 24,162,445
RAC: 0
20 million credit badge10 year member badge
Message 5425 - Posted: 9 Oct 2008, 20:51:58 UTC - in response to Message 5423.  

The qgaus call isn't in a loop, but the function that it is a part of could be. Chasing up, I'm finding a set of loops with a struct involved. I'm guessing around without the benefit of being able to debug and step through...

Yeah, it is not comfortable, but that's the route I also took.

Also, I do notice a difference between what I posted for gaussLegendre (gauleg) and what the project has after the while - the array positions. Again, without the capability to step it through, I don't know if it matters, but it appears to matter.

From a short look it appears one array starts at index "1", in the other example at index "0".
Edit: BTW, you are stating in bold two different figures, one of "180 million times", then another of "8 billion times". Without the benefit of being able to run the code in a debugger or have a profiler profile it, I can't definitively state that I think you're exaggerating things for dramatic effect, but clearly an approximate 44 times increase in the number of iterations that are deemed "excessive", particularly at that high of a number, seems unlikely. If it really goes through there "more than it should", then expressing that as either "more than it should" or with a number closer to what it really does would make your claim more believable... Over-dramatizing isn't helpful. Now if it really does loop 8 billion times, then yeah, there's a severe problem.

I am not exaggerating. These are correct numbers. It runs ~180 million times through qgaus (only counting the calculate_integrals part as it is 99.x% of the total execution time). The factor between the number of calls to qgaus and lbr2xyz is not 44 but 45, so it is definately more than 8 billion times (for a 260cr WU).
The number (for lbr2xyz) is approximately mu_steps * nu_steps * r_steps * convolve * 3
For the 260cr WUs (just look in the parameters file) the actual values are:
N (lbr2xyz) = 1600 * 80 * 700 * 30 * 3 = 8,064,000,000
These calls alone take longer than the whole WU with my app.

So you agree, there is a severe problem?
ID: 5425 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5426 - Posted: 9 Oct 2008, 22:24:30 UTC - in response to Message 5425.  


So you agree, there is a severe problem?


I don't know enough about it to say for sure. I've only taken a quick look. Where I worked before, we had to chase things all over creation sometimes, so I'm kinda used to that, but I've never been very good at nested loops with the question of "what is the output"? Maybe I have test-taking anxiety? Not sure...

I'll look at it some over the weekend. I assume that you have watched the variables at that point and they are the same as the original app through the whole process and/or the output is the same for both your app and the project's app?
ID: 5426 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5427 - Posted: 9 Oct 2008, 22:27:03 UTC
Last modified: 9 Oct 2008, 22:29:27 UTC

Changed the thread title. Hope everyone likes the nice positive spin ;-)

Also, since we're going into code discussion, the thread might be moved over there. I'm not going to request it myself, but if the mods request it or someone else requests it, I wouldn't be opposed...
ID: 5427 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Cluster Physik

Send message
Joined: 26 Jul 08
Posts: 627
Credit: 94,940,203
RAC: 0
50 million credit badge10 year member badgeextraordinary contributions badge
Message 5428 - Posted: 9 Oct 2008, 23:09:32 UTC - in response to Message 5426.  

I'll look at it some over the weekend. I assume that you have watched the variables at that point and they are the same as the original app through the whole process and/or the output is the same for both your app and the project's app?

You can check what values are passed to the functions. If they are the same, I would assume, the output is also the same every time ;)
And of course I've only cut out those calls which passed the same values as preceding ones. No magic here.
ID: 5428 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
vraa

Send message
Joined: 6 Mar 08
Posts: 8
Credit: 6,870,220
RAC: 0
5 million credit badge10 year member badge
Message 5430 - Posted: 10 Oct 2008, 1:56:11 UTC
Last modified: 10 Oct 2008, 1:57:18 UTC

http://transmission.pastebin.com/f69401f4

Sorry for the length
I was able to figure out how to run the Astronomy process through Shark, a popular profiling tool on Mac OS X
I am not a programmer though, so I'm not too familar on how to interpret the results
I hope the above is helpful in determining where the eyeballs should be focusing
I can produce many more if needed, this was a 30 second sample with every sample done every 1 ms.
ID: 5430 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
Brian Silvers

Send message
Joined: 21 Aug 08
Posts: 625
Credit: 558,425
RAC: 0
500 thousand credit badge10 year member badge
Message 5431 - Posted: 10 Oct 2008, 2:33:58 UTC - in response to Message 5430.  
Last modified: 10 Oct 2008, 2:43:14 UTC

I am not a programmer though, so I'm not too familar on how to interpret the results


I am a programmer, but not one that has ever really had to deal with performance profiling and tuning of an application. I see that it spent a lot of time in what I had gone back up to earlier today, calculate_integrals. It has pointers to what I believe are 3 structs being sent into it. I may be wrong about structs. It has been quite a while since actively dealing with C or C++. This is where I start disliking C, with all the pointers and referencing and dereferencing and manual memory management rather than the "easy" way of garbage collection, etc, etc, etc... However, I know that opinion comes from the fact that I do "business logic" applications, not scientific or math applications...

Anyway, if I'm feeling really industrious, sometime over the weekend I may try to get my old Visual Studio 6 installed...and see if I can remember how in the heck to compile / make the thing... .Net makes it real easy to just load the .prj or .sln...
ID: 5431 · Rating: 0 · rate: Rate + / Rate - Report as offensive     Reply Quote
1 · 2 · Next

Message boards : Number crunching : Source code improvement discussion

©2020 Astroinformatics Group