www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - std.range should support recursion (Was: One-line FFT, nice!)

reply "Mehrdad" <wfunction hotmail.com> writes:
I thought I'd take the opportunity to point this out:

The one-line FFT in D is pretty inefficient because it allocates 
memory.

If std.range supported recursion (i.e. by providing a different 
implementation for ranges that can be implemented without 
creating new times, i.e. Stride of Stride == Stride), then it 
would make the library a lot more usable and less bloated.


My one-line FFT illustrates it perfectly:


import std.algorithm;
import std.math;
import std.range;
typeof(R.init.stride(0)) dft(R)(R v)
{
	return v.length > 1 ?
	  (p => chain(map!(q => q[0] + q[1])(p), map!(q => q[0] - 
q[1])(p)))
	  (zip(dft(v.stride(2)),
	       map!(p => p[1] * expi(p[0] * -2 * PI / v.length))
	       (zip(iota(v.length / 2), dft(v.drop(1).stride(2)))))) : v;
}
void main() { dft([1.0, 2, 3]); }


Side note: the error messages are also hard to read.
Sep 25 2012
next sibling parent reply "Mehrdad" <wfunction hotmail.com> writes:
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
 without creating new times

new types**
Sep 25 2012
parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 9/25/12 4:23 AM, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
 without creating new times

new types**

Ah, better now. Still it would be great to explain it more :o). Andrei
Sep 25 2012
parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 9/25/12 11:42 AM, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei Alexandrescu wrote:
 On 9/25/12 4:23 AM, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
 without creating new times

new types**

Ah, better now. Still it would be great to explain it more :o). Andrei

Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3)))); otherwise recursion with these ranges is impossible.

I think all of the above are doable and useful. Please file a bug report containing these and any others you could reasonably think of. Thanks! Andrei
Sep 25 2012
prev sibling next sibling parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 9/25/12 4:22 AM, Mehrdad wrote:
 If std.range supported recursion (i.e. by providing a different
 implementation for ranges that can be implemented without creating new
 times, i.e. Stride of Stride == Stride), then it would make the library
 a lot more usable and less bloated.

I'm not sure I understand this, and it seems I should. Could you please explain (perhaps with a simpler example than FFT)? Thanks, Andrei
Sep 25 2012
prev sibling next sibling parent "Mehrdad" <wfunction hotmail.com> writes:
On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei 
Alexandrescu wrote:
 On 9/25/12 4:23 AM, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
 without creating new times

new types**

Ah, better now. Still it would be great to explain it more :o). Andrei

Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3)))); otherwise recursion with these ranges is impossible. The FFT example took the odd- and even-indexed numbers with stride(), but it couldn't recursively do this because the type system prevented it from doing so. So I was forced to copy the array unnecessarily every time. Also, foo should be implicitly convertible to typeof(foo.stride(1)), which also makes recursion easier.
Sep 25 2012
prev sibling next sibling parent "monarch_dodra" <monarchdodra gmail.com> writes:
On Tuesday, 25 September 2012 at 15:41:42 UTC, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei 
 Alexandrescu wrote:
 On 9/25/12 4:23 AM, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
 without creating new times

new types**

Ah, better now. Still it would be great to explain it more :o). Andrei

Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3))));

I can't comment on the rest of your points, but stride and take DO check for type recursivity, and drop always returns the same type as input anyways. Failure of ANY of these asserts is a bug. What where your inputs? //----------------------- import std.range; struct S { enum empty = false; void popFront(){}; property int front(){return 1;} } void main() { S foo; static assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); static assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); static assert(is(typeof(foo) == typeof(foo.drop(2)))); //Or this static assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3)))); } //-----------------------
Sep 25 2012
prev sibling next sibling parent "Mehrdad" <wfunction hotmail.com> writes:
On Tuesday, 25 September 2012 at 16:03:22 UTC, monarch_dodra 
wrote:
 On Tuesday, 25 September 2012 at 15:41:42 UTC, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei 
 Alexandrescu wrote:
 On 9/25/12 4:23 AM, Mehrdad wrote:
 On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
 without creating new times

new types**

Ah, better now. Still it would be great to explain it more :o). Andrei

Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3))));

I can't comment on the rest of your points, but stride and take DO check for type recursivity, and drop always returns the same type as input anyways. Failure of ANY of these asserts is a bug. What where your inputs?

I just wrote down the assert's on the fly, actually. Maybe I'm just misinterpreting the error then, and the problem is somewhere else? The code I was trying to compile is this (sorry it's ugly): import std.algorithm, std.math, std.range; typeof(R.init.stride(0)) dft(R)(R v) { return v.length <= 1 ? v.stride(1) : (p => chain(map!(q => q[0] + q[1])(p), map!(q => q[0] - q[1])(p))) (zip(dft(v.stride(2)), map!(p => p[1] * expi(p[0] * -2 * PI / v.length)) (zip(iota(v.length / 2), dft(v.drop(1).stride(2)))))); } void main() { dft([1.0, 2, 3]); } Which gives the following error: Test.d(5): Error: incompatible types for ((stride(v,1u)) ? ((*delegate system Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p) { return chain(map(p),map(p)); } )(zip(dft(stride(v,2u)),map(zip(iota(v.length() / 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result' Test.d(10): Error: template instance Test.dft!(Result) error instantiating Test.d:19: instantiated from here: dft!(double[]) Test.d(5): Error: incompatible types for ((stride(v,1u)) ? ((*delegate system Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p) { return chain(map(p),map(p)); } )(zip(dft(stride(v,2u)),map(zip(iota(v.length / 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result' Test.d(19): Error: template instance Test.dft!(double[]) error instantiating How should I interpret it? Thanks!
Sep 25 2012
prev sibling next sibling parent "Mehrdad" <wfunction hotmail.com> writes:
On Tuesday, 25 September 2012 at 16:31:15 UTC, Andrei 
Alexandrescu wrote:
 I think all of the above are doable and useful. Please file a 
 bug report containing these and any others you could reasonably 
 think of.

 Thanks!

 Andrei

Sure! I'll file them as soon as I know they're bugs -- monarch_dodra's comment makes me think it might just be my misinterpretation of the error.
Sep 25 2012
prev sibling next sibling parent "jerro" <a a.com> writes:
 import std.algorithm, std.math, std.range;
 typeof(R.init.stride(0)) dft(R)(R v)
 {
 	return v.length <= 1
 	? v.stride(1)
 	: (p =>
 	   chain(map!(q => q[0] + q[1])(p),
 	         map!(q => q[0] - q[1])(p)))
 	  (zip(dft(v.stride(2)),
 		map!(p => p[1] * expi(p[0] * -2 * PI / v.length))
 		(zip(iota(v.length / 2),
 		 dft(v.drop(1).stride(2))))));
 }
 void main() { dft([1.0, 2, 3]); }


 Which gives the following error:


 Test.d(5): Error: incompatible types for ((stride(v,1u)) ? 
 ((*delegate  system 
 Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p)
 {
 	return chain(map(p),map(p));
 }
 )(zip(dft(stride(v,2u)),map(zip(iota(v.length() / 
 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
 Test.d(10): Error: template instance Test.dft!(Result) error 
 instantiating
 Test.d:19: instantiated from here: dft!(double[])
 Test.d(5): Error: incompatible types for ((stride(v,1u)) ? 
 ((*delegate  system 
 Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p)
 {
 	return chain(map(p),map(p));
 }
 )(zip(dft(stride(v,2u)),map(zip(iota(v.length / 
 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
 Test.d(19): Error: template instance Test.dft!(double[]) error 
 instantiating



 How should I interpret it?

 Thanks!

One possible reason for this error could be that you are returning a result of chain() if length is larger than 1 and a result of stride() otherwise.
Sep 25 2012
prev sibling next sibling parent "jerro" <a a.com> writes:
On Tuesday, 25 September 2012 at 17:48:49 UTC, jerro wrote:
 import std.algorithm, std.math, std.range;
 typeof(R.init.stride(0)) dft(R)(R v)
 {
 	return v.length <= 1
 	? v.stride(1)
 	: (p =>
 	   chain(map!(q => q[0] + q[1])(p),
 	         map!(q => q[0] - q[1])(p)))
 	  (zip(dft(v.stride(2)),
 		map!(p => p[1] * expi(p[0] * -2 * PI / v.length))
 		(zip(iota(v.length / 2),
 		 dft(v.drop(1).stride(2))))));
 }
 void main() { dft([1.0, 2, 3]); }


 Which gives the following error:


 Test.d(5): Error: incompatible types for ((stride(v,1u)) ? 
 ((*delegate  system 
 Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) 
 p)
 {
 	return chain(map(p),map(p));
 }
 )(zip(dft(stride(v,2u)),map(zip(iota(v.length() / 
 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
 Test.d(10): Error: template instance Test.dft!(Result) error 
 instantiating
 Test.d:19: instantiated from here: dft!(double[])
 Test.d(5): Error: incompatible types for ((stride(v,1u)) ? 
 ((*delegate  system 
 Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) 
 p)
 {
 	return chain(map(p),map(p));
 }
 )(zip(dft(stride(v,2u)),map(zip(iota(v.length / 
 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
 Test.d(19): Error: template instance Test.dft!(double[]) error 
 instantiating



 How should I interpret it?

 Thanks!

One possible reason for this error could be that you are returning a result of chain() if length is larger than 1 and a result of stride() otherwise.

I'd like to add that I don't think you can make this work the way you meant it to. The problem is that you return a chain when the length is 2, a chain of chains when the length is 4, and so on. What fft of length n would actually need to return is a binary tree of ranges with n leaves. The size of memory needed for the value returned from fft therefore depends on the length of the range it was given as an argument. So you need to either use heap allocated memory for the return type, or the return type needs to depend on the parameter range's length, which would mean that the parameter range's length needs to be a template parameter. I think there is one thing in this code that will hurt performance much, much, more than allocations. This code will compute elements of the result lazily. So each time you want to read an element from the resulting range, O(log(n)) functions passed to map() will need to be computed. The problem is that each of those functions computes sine and cosine, so sine and cosine need to be computed O(log(n)) times for each element. To get all n elements, you will need to compute them O(n log(n)). Because computing sine and cosine is about two orders of magnitude slower than multiplication and division, this will be very slow.
Sep 25 2012
prev sibling next sibling parent "jerro" <a a.com> writes:
 I think there is one thing in this code that will hurt 
 performance much, much, more than allocations. This code will 
 compute elements of the result lazily. So each time you want to 
 read an element from the resulting range, O(log(n)) functions 
 passed to map() will need to be computed. The problem is that 
 each of those functions computes sine and cosine, so sine and 
 cosine need to be computed O(log(n)) times for each element. To 
 get all n elements, you will need to compute them O(n log(n)). 
 Because computing sine and cosine is about two orders of 
 magnitude slower than multiplication and division, this will be 
 very slow.

I was wrong about the complexity. Because each element of the result depends on all the elements of the argument range, you actually need O(n) function calls to compute each element of the result and O(n*n) function calls(and sine and cosine computations) to compute all of them. You would need to use memoization to get reasonable complexity.
Sep 25 2012
prev sibling parent "Mehrdad" <wfunction hotmail.com> writes:
On Tuesday, 25 September 2012 at 18:33:45 UTC, jerro wrote:
 I think there is one thing in this code that will hurt 
 performance much, much, more than allocations. This code will 
 compute elements of the result lazily. So each time you want 
 to read an element from the resulting range, O(log(n)) 
 functions passed to map() will need to be computed. The 
 problem is that each of those functions computes sine and 
 cosine, so sine and cosine need to be computed O(log(n)) times 
 for each element. To get all n elements, you will need to 
 compute them O(n log(n)). Because computing sine and cosine is 
 about two orders of magnitude slower than multiplication and 
 division, this will be very slow.

I was wrong about the complexity. Because each element of the result depends on all the elements of the argument range, you actually need O(n) function calls to compute each element of the result and O(n*n) function calls(and sine and cosine computations) to compute all of them. You would need to use memoization to get reasonable complexity.

Great point, I hadn't really thought about the laziness before -- I was doing this in Python and just thought it might be fun to translate it to D. :P
 I'd like to add that I don't think you can make this work the 
 way you meant it to. The problem is that you return a chain 
 when the length is 2, a chain of chains when the length is 4, 
 and so on. What fft of length n would actually need to return 
 is a binary tree of ranges with n leaves. The size of memory 
 needed for the value returned from fft therefore depends on the 
 length of the range it was given as an argument. So you need to 
 either use heap allocated memory for the return type, or the 
 return type needs to depend on the parameter range's length, 
 which would mean that the parameter range's length needs to be 
 a template parameter.

Ohh huh... sounds like you're right, lemme think about it a bit more though. Thanks! :)
Sep 25 2012