Sign in to follow this  
jmau0438

Trouble with sqrt & template functions

Recommended Posts

Hello out there, I'm writting a stats utility header filled with templated functions, and I've run into some trouble when I try using the sqrt function in <math.h>. Here is the function:
template < typename TYPE > TYPE STD_DEVIATION ( const TYPE *Array, const int Length ) {

		//Declare & Initialize the mean to the average of the sample
		TYPE Mean = MEAN < TYPE > ( Array, Length );

		//Declare an array containing ( X - M ) data
		TYPE *X_M = new TYPE [ Length ];

		//Declare an array containing ( X - M ) 2 data
		TYPE *X_M_SQRT = new TYPE [ Length ];

		//Declare a variable to store the sum of ( X - m ) 2 & Initialize to zero
		TYPE X_M_SQRT_SUM = 0;

		//For each value within a sample, calculate & store the mean & deviation by...
		for ( int intLoop = 0; intLoop < Length; intLoop++ ) {

			//...subtracting the mean from the current sample value & store it in the ( X - M ) array...
			X_M [ intLoop ] = ( Array [ intLoop ] - Mean );

			//...and squaring the result of the previous calculation & store it in the ( X - M ) 2 array...
			X_M_SQRT [ intLoop ] = ( sqrt ( X_M [ intLoop ] ) );

		};

		//Calculate the sum of ( X - M ) 2
		for ( int intLoop = 0; intLoop < Length; intLoop++ ) { X_M_SQRT_SUM += X_M_SQRT [ intLoop ]; };

		//Purge the ( X - M ) & ( X - M ) 2 arrays from memory
		delete X_M, X_M_SQRT;
		X_M, X_M_SQRT = 0;

		//Return the standard deviation
		return ( sqrt ( X_M_SQRT_SUM ) / sqrt ( Length - 1 ) );

	};


HERE IS THE ERROR error C2668: 'sqrt' : ambiguous call to overloaded function I've tried hacking in some type casting tricks, but i usually end up with a -1.#IND result, which I'm pretty sure is a type of NAN value. This would make sense if I was trying to square zero or maybe even a negative number, but I know I'm not. I know much of the problem revolves around the fact that I'm using a templated funtion and the compiler doesn't know the type during compilation, which is giving me this error. Is what I'm trying just not possible, or is there a work around? [Edited by - Zahlman on August 19, 2009 9:10:36 AM]

Share this post


Link to post
Share on other sites
The compiler knows the type, because templates aren't actually "compiled" when the compiler sees them. Rather, the compiler just holds onto the code, and it's only when the template is "instantiated" that the compiler fills in all the blanks and can generate code.

The problem is that you're supplying an integer to the sqrt function (when you perform sqrt(Length - 1)). Now, sqrt can take a float or a double. The compiler needs to convert your int into a float or a double - but which one? The compiler can't make this decision, and so it gives you an "ambiguous call" error.

Share this post


Link to post
Share on other sites
Never mind, I just found out that I was squaring a negative number. Oops!! Anyway, thanks so much for your help Sc4Freak, I really appreciate it!!

Share this post


Link to post
Share on other sites
Quote:
Original post by jmau0438
Never mind, I just found out that I was squaring a negative number. Oops!! Anyway, thanks so much for your help Sc4Freak, I really appreciate it!!

Yes, on this line:

X_M_SQRT [ intLoop ] = ( sqrt ( X_M [ intLoop ] ) );

You're supposed to square the number, not square root.

There are some other things I'd like to mention as a sidenote, however:
1) Use the [source][/source] tags to format your code nicely
2) You're doing a lot of unneccessary calculation. The following code does the same thing and is a lot more clearer and compact:
template<typename T>
T StdDev(const T* data, size_t num)
{
const T* end = data + num;

// If you don't know what std::accumulate does, it just adds stuff in
// the given range
T mean = std::accumulate(data, end, 0) / num;

T sum = 0;
for(; data != end; ++data)
{
T diff = *data - mean;
sum += diff * diff;
}
sum /= num;

return static_cast<T>(sqrt(static_cast<double>(sum)));
}

No memory allocation needed, and just one loop.
3) You leak memory. You allocate the arrays X_M and X_M_SQRT, but you call delete. You need to call delete[] if you allocated using new[].

Share this post


Link to post
Share on other sites
Quote:
Original post by jmau0438
//Purge the ( X - M ) & ( X - M ) 2 arrays from memory
delete X_M, X_M_SQRT;
X_M, X_M_SQRT = 0;


Unfortunately, these lines do not do what you expect. They're equivalent to the following:


delete X_M;

X_M_SQRT; // Does nothing

X_M; // Does nothing

X_M_SQRT = 0;




You have to perform each action individually, as follows (incorporating Sc4Freak's delete[] fix):


delete[] X_M, delete[] X_M_SQRT;
X_M = 0, X_M_SQRT = 0;




The comma separates statements (just like the semi-colon). It does not apply one operator to multiple values.

Share this post


Link to post
Share on other sites
Hi

Also, you don't in general need to put semicolons after a set of curly braces - only if it's a struct/class definition.



for(int i=0;i<maxloop;++i)
{
//Do stuff
//Do more stuff
}



will work fine. Curly brackets are meant to group statements into a block - they aren't statements in themselves, so they don't need to be finished with the semicolon.

Share this post


Link to post
Share on other sites
Quote:
delete X_M, X_M_SQRT;
As said, this doesn't work.

Unless you use the improved version with no allocation, use something like this instead:
//Declare an array containing ( X - M ) data
std::vector<TYPE> X_M(length);
//Declare an array containing ( X - M ) 2 data
std::vector<TYPE> X_M_SQRT(length);

Share this post


Link to post
Share on other sites
Quote:
Original post by Sc4Freak
3) You leak memory. You allocate the arrays X_M and X_M_SQRT, but you call delete. You need to call delete[] if you allocated using new[].


That's not a memory leak; it's undefined behaviour.

And as long as we're showing off templates and the standard library algorithms, why not adapt the interface to follow the standard library conventions as well? :)

Oh, and in case the input type is an integer, you might not want to be doing integer division all over. Fortunately, there is an algebraic simplication which gets around that...

When we calculate the squared sum, we calculate the sum of all (x - x0)2, i.e. the sum of all (x2 - 2x*x0 + x02). Those sums are separable, and the sum of 2x*x0 is 2x0 times the sum of the x values; but since x0 is the mean value, the sum of the x values is simply x0 * N.

So we end up with the (sum of x2) - 2N x02 + N x02, or just (sum of x2) - N x02. And since (sum of x) == (N x0), again, we can simplify as (sum of x2) - ((sum of x)2 / N). :)

Finally, it doesn't really make sense to force the result type back to the same as the element type. The standard deviation is intrinsically a decimal quantity.


template <typename ForwardIterator> double std_deviation(ForwardIterator begin, ForwardIterator end) {
double N = static_cast<double>(std::distance(begin, end));

typedef typename ForwardIterator::value_type T;
T sum = std::accumulate(begin, end, 0);

T sum_of_squares = 0;
for (; begin != end; ++begin) {
sum_of_squares += *begin * *begin;
}

// Now all the division comes at the end.
return sqrt(
static_cast<double>(sum_of_squares) -
static_cast<double>(sum * sum) / N
) / sqrt(N - 1);
}

Share this post


Link to post
Share on other sites
Wow, you guys are awesome. I guess I respond in parts:

Sc4Freak - Square not Square Root

I was to focused on the square root problem. After you showed me how to fix the sqrt() function call, it took but a few seconds for me to realize that when I performed the square root operation (which is flat wrong at this stage in the formula) that I was creating negative numbers, which gave me the NAN value (1.#IND).

Sc4Freak, scjohnno, Zahlman, & Antheus - Memory Leak

Ya, I caught the memory leak too. Again, I was way to focused on the sqrt() problem so I didn't catch it until last night. Sloppy I know, but hey, it was 4 in the morning. I simplified the code a bit so I'm only using a single local array, and replaced the previous deletion with this:

//Purge the ( X - M ) array from memory
delete [] X_M, X_M = 0;

shaolinspin - Semi-colon after curly brace

I know about the semi-colon rule, but it was a bad habit I picked up in college. To this day I habitually throw a semi-colon behind closing curly braces. It haven't experienced any trouble caused by the habit, but just out of curiosity are there any adverse effects to adding the semi-colon?

Share this post


Link to post
Share on other sites
I don't think so, but someone more knowledgeable than I might beg to differ. As far as I'm aware, an extra semi-colon will result in an empty statement in the assembled code, that the compiler will probably optimise away for you.

Share this post


Link to post
Share on other sites
Important tip: the comma operator is a subtle beast. Don't try to be clever, you'll only get burned.

Just write it on separate lines. Or use std::vector<>. Why aren't you using std::vector<>?

delete [] X_M;
X_M = 0;

Share this post


Link to post
Share on other sites
I actually am using vectors. What I had posted was kind of like a draft. I took in almost all the advise of the fine programmers at game dev. For those who are curious of my final solution, I suppose I'll post it here:

template < typename TYPE > TYPE STD_DEVIATION ( const vector < TYPE > Array ) {

//The sum of all values in the sample/array, or the sum of X
TYPE SampleSum = 0;
//The squared sum of all values in the sample/array, or the sum of (X)2
TYPE SampleSumSquared = 0;
//Constant value representing the number of (actively used) elements in the vector.
const int Length = Array.size();

//Compute & store the sum of the samples' raw & squared values
for ( int intLoop = 0; intLoop < Length; intLoop++ ) {
SampleSum += Array [ intLoop ];
SampleSumSquared += ( Array [ intLoop ] * Array [ intLoop ] );
};

//Return the population standard deviation
return sqrt ( ( SampleSumSquared - ( SampleSum * SampleSum ) /
Length ) / ( Length - 1 ) );
};

What do you fellas think?

Share this post


Link to post
Share on other sites
Looks good. And yeah, as long as you need the for-loop anyway for the sum of squares (it could be done with std::accumulate, std::transform and a couple other pieces, I think, but it really doesn't do a very good job of looking any simpler :) ), might as well use it for the sum of values too. ;)

Only things I would change:

- Pass the vector in by reference, and don't call it 'Array' when it isn't one. :) (Actually, instead of trying to describe what kind of container it is, better to say a little about what its contents are. ;) )

- Iterate using an iterator, and save the value in a temporary:


// BTW, a loop iteration variable doesn't benefit from that 'int' prefix
// any more than anything else. :)
typedef vector<TYPE>::const_iterator const_iterator;
for (const_iterator it = data.begin(); it != data.end(); ++it) {
const TYPE& value = *it;
SampleSum += value;
SampleSumSquared += (value * value);
} // no semicolon needed here


- It really doesn't need that much commenting ;)

Share this post


Link to post
Share on other sites
I used accumulate in the function, and I gotta say it made it way more readable, check it out:

template < typename TYPE > TYPE STD_DEVIATION ( const vector < TYPE > &Vec ) {

//The sum of all values in the sample/array, or the sum of X
TYPE SampleSum = accumulate ( Vec.begin(), Vec.end(), 0 );
//The squared sum of all values in the sample/array, or the sum of (X)2
TYPE SampleSumSquared = accumulate ( Vec.begin(), Vec.end(), 0, AccumulateSquaredSum < TYPE > );
//Constant value representing the number of (actively used) elements in the vector.
const int Length = Vec.size();

//Return the population standard deviation
return sqrt ( ( SampleSumSquared - ( SampleSum * SampleSum ) / Length ) / ( Length - 1 ) );

};

Yea I kept the comments, but I can't tell you how many times I wrote something nice and just used it for like a year, and when I wanted to modify it I ended up spending more time than I would have liked understanding what it was that I did. At least this way I won't forget, and explaining it helps me to better understand what it was that I did. I know we've probably beat this whole thing to death, but what do ya think?

Share this post


Link to post
Share on other sites
I think you should drop the all-caps, they are hard to read and they are universally understood to stand for constant values, not types or function names.

I generally wouldn't mark a variable like "Length" const. While technically correct it can add unnecessary mental processing because most programmers associate "const" with pointers or reference types (and constant values of course). This is a minor point, a more important one is consistency. If Length is constant - why not SampleSum or SampleSumSquared?

Quote:

Yea I kept the comments, but I can't tell you how many times I wrote something nice and just used it for like a year, and when I wanted to modify it I ended up spending more time than I would have liked understanding what it was that I did. At least this way I won't forget, and explaining it helps me to better understand what it was that I did. I know we've probably beat this whole thing to death, but what do ya think?

The *intent* behind your is fine, but they are a little verbose. The don't actually give any more information than the variable names (which are well chosen here) and algorithms do. This is what we call self documenting code, to anyone familiar with the standard algorithms the code is clear.

I would me more concerned with documenting the rather obscure looking formula used in the return value than the variables being input into it. Even including a link you can refer to later would be a good idea.

Also, I would be tempted to break it over a few lines, like Zahlman did, or introduce a few temporaries. Large expressions are difficult for the mind to parse - make it easier.

Share this post


Link to post
Share on other sites
Quote:
Original post by Zahlman
So we end up with the (sum of x2) - 2N x02 + N x02, or just (sum of x2) - N x02. And since (sum of x) == (N x0), again, we can simplify as (sum of x2) - ((sum of x)2 / N). :)


Unfortunately, floating point arithmetic doesn't follow all the same laws as ideal mathematics. This transformation is quite literally a textbook example of a numerically unstable algorithm (ex: pg 19 of "Scientific Computing An Introduction Survey" by Heath). Taking the difference of two sums that can be very close to one another can result in what is called catastrophic cancellation and can cause standard deviation calculations orders of magnitude different than the actual value and, on occasion, taking the square root of a negative number.

Wikipedia lists a few more stable algorithms for calculating variance.

Share this post


Link to post
Share on other sites
And here I was doing it specifically to avoid round-off error by delaying the division when the template type is integral... :P

Anyway, as a rule of thumb, if it isn't readable without the comments, it isn't really readable with them, either. Comments aren't there to enhance readability so much as to provide direction to the person reading the code.

That doesn't always mean you have to rewrite it again. It sometimes means you have to trust your ability to understand (and write) code a little more. :)

Share this post


Link to post
Share on other sites
Thanks for your posting guys. Again, I guess I'll have to answer in parts:

Zahlman & Rip-Off - Comments

In terms of commenting, I'm not trying to improve readability so much as to leave behind bread-crumbs for later. As Zahlman said, "Comments aren't there to enhance readability so much as to provide direction to the person reading the code." This has kind of been my aim with comments, to leave behind direction. Maybe I'm not hitting that up as well as I could, so I'll work on it. Again to quote Zahlman, "It sometimes means you have to trust your ability to understand (and write) code a little more.", and maybe trusting myself in this respect is what I'm really having a problem with.

Rip-Off - Constant Consistency and the Return Value

I agree on your point of consistency, so the "SampleSum" & "SampleSumSquared" are now constant. Although I understand what you mean in how a constant cast can "add unnecessary mental processing because most programmers associate "const" with pointers or reference types (and constant values of course)", I've chosen to retain the constant cast as identifying them as simply constants isn't too far of a strech in this case. Your point on further commenting on the return value is also a good one, and I'll make changes to comment after this post.

SiCrane - Statistics & Negative Numbers

What you say in regards to negative numbers in this formula is absolutely true. However, I have found in most (keyword is most here) statistical calculations when you get a negative result you should check your math before anything else. In fact, in this very discussion thread I had a problem were I was performing a squared root of a value when I should have just squared it, which yielded a negative number just prior to the final algorithm. In this case, I got C/C++'s NaN value (-1.#IND), which did indeed suck. However, it was my math, not the algorithm, which got the better of me. Your point is still valid, but in this case (that being statistical math) I'd prefer it broke here as it will indicate to me that my math is wrong somewhere. Man, I should write an exception for just that! Awesome!

Share this post


Link to post
Share on other sites
Quote:
Original post by jmau0438
Your point is still valid, but in this case (that being statistical math) I'd prefer it broke here as it will indicate to me that my math is wrong somewhere.

That makes no sense. If it breaks here it doesn't indicate anything about your math anywhere else. It only indicates that you used a numerically unstable algorithm here.

Share this post


Link to post
Share on other sites
SiCrane -

I followed the link you left me to wikipedia. It covered several alternatives, but most had secondary risks. There was one that caught my eye, so could you tell me it this is it?

Taken from wikipedia.org

[Edited by - jmau0438 on August 27, 2009 10:36:13 AM]

Share this post


Link to post
Share on other sites
Trying to post an image with the source being your hard drive isn't something that's going to work happily.

Share this post


Link to post
Share on other sites
Yea I found that out, sorry man I'm new to this. It's an image of the formula above the Weighted incremental algorithm header.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this