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 reply "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
next sibling parent reply "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
parent reply "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
parent reply "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
parent reply "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
parent reply "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
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
prev sibling parent reply 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
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 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