comp.lang.idl-pvwave archive - RDF feed
type: channel
description: Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
items: Accelerating a one-line program doing matrix multiplication
items: Re: Accelerating a one-line program doing matrix multiplication

Accelerating a one-line program doing matrix multiplication
type: item
description: Hi all, I wrote a one-line function to convert a list of points from "voxel coordinates" (image coordinates) to "real coordinates" (physical coordinates): ;input: the points "vc", the spatial origin of an image v0 and its x, y, and z orientation vectors (v1,v2,v3). FUNCTION vc2rc, v0,v1,v2,v3,vc RETURN, [[v1],[v2],[v3]] # vc + REBIN(v0, SIZE(vc, /DIMENSIONS)) END For example, I give the image coordinate [8,1,0] and I want as output something like [34.25, 4.12, 0], indicating the location of this voxel in space. And the same thing but, instead of having one input point, having several millions. The function looks simple to me and it works great. BUT, for large images (e.g. 500x500x200 voxels), it is terribly slow and uses way too much memory... Am I doing something wrong, could I save speed somewhere? I guess there should be some way to accelerate it, but I am not able to see how... I also have the opposite function, in my opinion also too slow (though faster than the other)... FUNCTION rc2vc_round, v0,v1,v2,v3,rc RETURN, ROUND((rc - REBIN(v0, SIZE(rc, /DIMENSIONS))) ## INVERT([[v1], [v2],[v3]])) END I would be really grateful for any clue!
creator: Axel Martínez
date: 2010-09-27T09:18:43-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 29, 4:49 pm, chris <rog...@googlemail.com> wrote: > On 29 Sep., 12:48, Axel M <axe...@gmail.com> wrote: > >> I have to admit that I did not understand this proposed use of >> REPLICATE: > >>>> sometimes (replicate({temp:input},newsize)).(0) is faster then rebin > >> But it brought me a related question in mind: does IDL have a >> "REPLICATE" function for vectors instead of scalar values? I am using >> REBIN, but REBIN is thought for more advanced uses and probably >> suboptimal for a "replicate-like" use... right? > > Yes, it has as I mentioned uncommented above: > > sometimes (replicate({temp:input},newsize)).(0) is faster then rebin > > -> this means: > > IDL> a=findgen(3) > IDL> print,a > 0.000000 1.00000 2.00000 > IDL> print, rebin(a,3,5) > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > IDL> print, (replicate({temp:a},5)).(0) > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > > Regards > > CR Thanks!! Great, I did not know about this construction, and honestly I do not understand how it works (is there any documentation about it?). Anyways, I tried it, and unfortunately I saw that it needed ~20% longer (the complete function, not the rebin only). So, it is not faster.. but it is great though.
creator: Axel Martínez
date: 2010-09-29T15:24:23-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: > > [skip] > > sometimes (replicate({temp:input},newsize)).(0) is faster then rebin > > -> this means: > > IDL> a=findgen(3) > IDL> print,a > 0.000000 1.00000 2.00000 > IDL> print, rebin(a,3,5) > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > IDL> print, (replicate({temp:a},5)).(0) > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > 0.000000 1.00000 2.00000 > > Regards > > CR Another option is: a#(fltarr(5)+1) But I would expect this to be slower... Ciao, Paolo
creator: pgrigis
date: 2010-09-29T15:34:29-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > Great, I did not know about this construction, and honestly I do not > understand how it works (is there any documentation about it?). > Anyways, I tried it, and unfortunately I saw that it needed ~20% > longer (the complete function, not the rebin only). So, it is not > faster.. but it is great though. It is replicating a structure of a single field which contains the array input ({temp:input}), then selecting only a single field (the first, 0) of the resulting structure array. Documentation for this would be on creation and use of structures.
creator: penteado
date: 2010-09-29T15:45:56-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On 29 Sep., 17:45, Paulo Penteado <pp.pente...@gmail.com> wrote: > On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > >> Great, I did not know about this construction, and honestly I do not >> understand how it works (is there any documentation about it?). >> Anyways, I tried it, and unfortunately I saw that it needed ~20% >> longer (the complete function, not the rebin only). So, it is not >> faster.. but it is great though. > > It is replicating a structure of a single field which contains the > array input ({temp:input}), then selecting only a single field (the > first, 0) of the resulting structure array. Documentation for this > would be on creation and use of structures. Ok, I got it. Thanks! Then probably it is the memory allocation for the array of structures which takes so long... it would be great if the ITT people would develop a _fast_ vector replicate, I fear rebinning is not the best option. In any case, based on the answers, I assume that my problem is rather on the matrix multiplication part, so I can probably do nothing for that. Thanks a lot
creator: Axel Martínez
date: 2010-09-29T15:55:33-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 29, 11:55 am, Axel M <axe...@gmail.com> wrote: > On 29 Sep., 17:45, Paulo Penteado <pp.pente...@gmail.com> wrote: > >> On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > >>> Great, I did not know about this construction, and honestly I do not >>> understand how it works (is there any documentation about it?). >>> Anyways, I tried it, and unfortunately I saw that it needed ~20% >>> longer (the complete function, not the rebin only). So, it is not >>> faster.. but it is great though. > >> It is replicating a structure of a single field which contains the >> array input ({temp:input}), then selecting only a single field (the >> first, 0) of the resulting structure array. Documentation for this >> would be on creation and use of structures. > > Ok, I got it. Thanks! Then probably it is the memory allocation for > the array of structures which takes so long... it would be great if > the ITT people would develop a _fast_ vector replicate, I fear > rebinning is not the best option. > > In any case, based on the answers, I assume that my problem is rather > on the matrix multiplication part, so I can probably do nothing for > that. > > Thanks a lot well considering your original problem - you need to apply a linear transformation to N vectors v_i=(x_i,y_i,z_i), for i going from 0 to a large N, right? I would just explicitely compute the transformed vectors z_i=(xx_i,yy_i,zz_i) by just writing out in the program the computation for every component, i.e. xx=x*c1+y*c2+z*c3+c4 and same for yy,zz with appropriate constant coefficients c1,c2,c3,c4 (that are the same for all i). But then maybe i misunderstood the problem... Ciao, Paolo
creator: pgrigis
date: 2010-09-29T16:05:58-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 29, 10:05 am, Paolo <pgri...@gmail.com> wrote: > On Sep 29, 11:55 am, Axel M <axe...@gmail.com> wrote: > > > >> On 29 Sep., 17:45, Paulo Penteado <pp.pente...@gmail.com> wrote: > >>> On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > >>>> Great, I did not know about this construction, and honestly I do not >>>> understand how it works (is there any documentation about it?). >>>> Anyways, I tried it, and unfortunately I saw that it needed ~20% >>>> longer (the complete function, not the rebin only). So, it is not >>>> faster.. but it is great though. > >>> It is replicating a structure of a single field which contains the >>> array input ({temp:input}), then selecting only a single field (the >>> first, 0) of the resulting structure array. Documentation for this >>> would be on creation and use of structures. > >> Ok, I got it. Thanks! Then probably it is the memory allocation for >> the array of structures which takes so long... it would be great if >> the ITT people would develop a _fast_ vector replicate, I fear >> rebinning is not the best option. > >> In any case, based on the answers, I assume that my problem is rather >> on the matrix multiplication part, so I can probably do nothing for >> that. > >> Thanks a lot > > well considering your original problem - you need to apply > a linear transformation to N vectors v_i=(x_i,y_i,z_i), > for i going from 0 to a large N, right? > > I would just explicitely compute the transformed vectors > > z_i=(xx_i,yy_i,zz_i) > > by just writing out in the program the computation for every > component, > i.e. > > xx=x*c1+y*c2+z*c3+c4 > and same for yy,zz with appropriate constant coefficients c1,c2,c3,c4 > (that are the same for all i). > > But then maybe i misunderstood the problem... > > Ciao, > Paolo Yeah, I think you are right. Another way to see it: FUNCTION vc2rc, v0,v1,v2,v3,vc xform = [[v1],[v2],[v3]] n = <number of points in vc> for i=0, n-1 temp = vc[*,i] temp = temp # xform + v0 vc[*,i] = temp end END Thi…
creator: Karl[1]
date: 2010-09-29T16:57:05-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 29, 6:57 pm, Karl <karl.w.schu...@gmail.com> wrote: > On Sep 29, 10:05 am, Paolo <pgri...@gmail.com> wrote: > > > >> On Sep 29, 11:55 am, Axel M <axe...@gmail.com> wrote: > >>> On 29 Sep., 17:45, Paulo Penteado <pp.pente...@gmail.com> wrote: > >>>> On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > >>>> > Great, I did not know about this construction, and honestly I do not >>>> > understand how it works (is there any documentation about it?). >>>> > Anyways, I tried it, and unfortunately I saw that it needed ~20% >>>> > longer (the complete function, not the rebin only). So, it is not >>>> > faster.. but it is great though. > >>>> It is replicating a structure of a single field which contains the >>>> array input ({temp:input}), then selecting only a single field (the >>>> first, 0) of the resulting structure array. Documentation for this >>>> would be on creation and use of structures. > >>> Ok, I got it. Thanks! Then probably it is the memory allocation for >>> the array of structures which takes so long... it would be great if >>> the ITT people would develop a _fast_ vector replicate, I fear >>> rebinning is not the best option. > >>> In any case, based on the answers, I assume that my problem is rather >>> on the matrix multiplication part, so I can probably do nothing for >>> that. > >>> Thanks a lot > >> well considering your original problem - you need to apply >> a linear transformation to N vectors v_i=(x_i,y_i,z_i), >> for i going from 0 to a large N, right? > >> I would just explicitely compute the transformed vectors > >> z_i=(xx_i,yy_i,zz_i) > >> by just writing out in the program the computation for every >> component, >> i.e. > >> xx=x*c1+y*c2+z*c3+c4 >> and same for yy,zz with appropriate constant coefficients c1,c2,c3,c4 >> (that are the same for all i). > >> But then maybe i misunderstood the problem... > >> Ciao, >> Paolo > > Yeah, I think you are right. > > Another way to see it: > > FUNCTION vc2rc, v0,v1,v2,v3,vc > xform = [[v1],[v2…
creator: Axel Martínez
date: 2010-09-30T08:39:59-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 30, 10:39 am, Axel M <axe...@gmail.com> wrote: > On Sep 29, 6:57 pm, Karl <karl.w.schu...@gmail.com> wrote: > > > >> On Sep 29, 10:05 am, Paolo <pgri...@gmail.com> wrote: > >>> On Sep 29, 11:55 am, Axel M <axe...@gmail.com> wrote: > >>>> On 29 Sep., 17:45, Paulo Penteado <pp.pente...@gmail.com> wrote: > >>>> > On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > >>>> > > Great, I did not know about this construction, and honestly I do not >>>> > > understand how it works (is there any documentation about it?). >>>> > > Anyways, I tried it, and unfortunately I saw that it needed ~20% >>>> > > longer (the complete function, not the rebin only). So, it is not >>>> > > faster.. but it is great though. > >>>> > It is replicating a structure of a single field which contains the >>>> > array input ({temp:input}), then selecting only a single field (the >>>> > first, 0) of the resulting structure array. Documentation for this >>>> > would be on creation and use of structures. > >>>> Ok, I got it. Thanks! Then probably it is the memory allocation for >>>> the array of structures which takes so long... it would be great if >>>> the ITT people would develop a _fast_ vector replicate, I fear >>>> rebinning is not the best option. > >>>> In any case, based on the answers, I assume that my problem is rather >>>> on the matrix multiplication part, so I can probably do nothing for >>>> that. > >>>> Thanks a lot > >>> well considering your original problem - you need to apply >>> a linear transformation to N vectors v_i=(x_i,y_i,z_i), >>> for i going from 0 to a large N, right? > >>> I would just explicitely compute the transformed vectors > >>> z_i=(xx_i,yy_i,zz_i) > >>> by just writing out in the program the computation for every >>> component, >>> i.e. > >>> xx=x*c1+y*c2+z*c3+c4 >>> and same for yy,zz with appropriate constant coefficients c1,c2,c3,c4 >>> (that are the same for all i). > >>> But then maybe i misunderstood the problem... > >>> Ciao, >>> Paolo > >> Yeah, …
creator: Axel Martínez
date: 2010-09-30T09:41:19-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 30, 11:41 am, Axel M <axe...@gmail.com> wrote: > On Sep 30, 10:39 am, Axel M <axe...@gmail.com> wrote: > > > >> On Sep 29, 6:57 pm, Karl <karl.w.schu...@gmail.com> wrote: > >>> On Sep 29, 10:05 am, Paolo <pgri...@gmail.com> wrote: > >>>> On Sep 29, 11:55 am, Axel M <axe...@gmail.com> wrote: > >>>> > On 29 Sep., 17:45, Paulo Penteado <pp.pente...@gmail.com> wrote: > >>>> > > On Sep 29, 12:24 pm, Axel M <axe...@gmail.com> wrote: > >>>> > > > Great, I did not know about this construction, and honestly I do not >>>> > > > understand how it works (is there any documentation about it?). >>>> > > > Anyways, I tried it, and unfortunately I saw that it needed ~20% >>>> > > > longer (the complete function, not the rebin only). So, it is not >>>> > > > faster.. but it is great though. > >>>> > > It is replicating a structure of a single field which contains the >>>> > > array input ({temp:input}), then selecting only a single field (the >>>> > > first, 0) of the resulting structure array. Documentation for this >>>> > > would be on creation and use of structures. > >>>> > Ok, I got it. Thanks! Then probably it is the memory allocation for >>>> > the array of structures which takes so long... it would be great if >>>> > the ITT people would develop a _fast_ vector replicate, I fear >>>> > rebinning is not the best option. > >>>> > In any case, based on the answers, I assume that my problem is rather >>>> > on the matrix multiplication part, so I can probably do nothing for >>>> > that. > >>>> > Thanks a lot > >>>> well considering your original problem - you need to apply >>>> a linear transformation to N vectors v_i=(x_i,y_i,z_i), >>>> for i going from 0 to a large N, right? > >>>> I would just explicitely compute the transformed vectors > >>>> z_i=(xx_i,yy_i,zz_i) > >>>> by just writing out in the program the computation for every >>>> component, >>>> i.e. > >>>> xx=x*c1+y*c2+z*c3+c4 >>>> and same for yy,zz with appropriate constant coefficients c1,c2,c3,c4 >>>> (that …
creator: Axel Martínez
date: 2010-09-30T11:05:36-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: [skip] > > FUNCTION vc2rc_accel, v0,v1,v2,v3,vc > npoints = (SIZE(vc, /DIMENSIONS))[1] > for i=0L, npoints-1 DO BEGIN > vc[*,i] = vc[0,i] * v1 + vc[1,i] * v2 + vc[2,i] * v3 + v0 > endfor > RETURN, vc > END No, that's using a for loop - that's why it is slow. You want something like this (no loops): vc0=vc[0,*] vc1=vc[1,*] vc2=vc[2,*] vc[0,*]=vc0*v1[0]+vc1*v2[0]+vc3*v3[0] vc[1,*]=vc0*v1[1]+vc1*v2[1]+... vc[2,*]=vc0*v1[2]+... Ciao, Paolo
creator: pgrigis
date: 2010-09-30T14:59:59-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 30, 10:59 am, Paolo <pgri...@gmail.com> wrote: > [skip] > > > >> FUNCTION vc2rc_accel, v0,v1,v2,v3,vc >> npoints = (SIZE(vc, /DIMENSIONS))[1] >> for i=0L, npoints-1 DO BEGIN >> vc[*,i] = vc[0,i] * v1 + vc[1,i] * v2 + vc[2,i] * v3 + v0 >> endfor >> RETURN, vc >> END > > No, that's using a for loop - that's why it is slow. > You want something like this (no loops): > > vc0=vc[0,*] > vc1=vc[1,*] > vc2=vc[2,*] > vc[0,*]=vc0*v1[0]+vc1*v2[0]+vc3*v3[0] > vc[1,*]=vc0*v1[1]+vc1*v2[1]+... > vc[2,*]=vc0*v1[2]+... > > Ciao, > Paolo well, there's a bit of an index mismatch in there... vc[0,*]=vc0*v1[0]+vc1*v2[0]+vc2*v3[0] vc[1,*]=vc0*v1[1]+vc1*v2[1]+vc2*v3[1] vc[2,*]=vc0*v1[2]+vc1*v2[2]+vc2*v3[2]
creator: pgrigis
date: 2010-09-30T15:03:43-00:00

Re: Accelerating a one-line program doing matrix multiplication
type: item
description: On Sep 30, 5:03 pm, Paolo <pgri...@gmail.com> wrote: > On Sep 30, 10:59 am, Paolo <pgri...@gmail.com> wrote: > > > >> [skip] > >>> FUNCTION vc2rc_accel, v0,v1,v2,v3,vc >>> npoints = (SIZE(vc, /DIMENSIONS))[1] >>> for i=0L, npoints-1 DO BEGIN >>> vc[*,i] = vc[0,i] * v1 + vc[1,i] * v2 + vc[2,i] * v3 + v0 >>> endfor >>> RETURN, vc >>> END > >> No, that's using a for loop - that's why it is slow. >> You want something like this (no loops): > >> vc0=vc[0,*] >> vc1=vc[1,*] >> vc2=vc[2,*] >> vc[0,*]=vc0*v1[0]+vc1*v2[0]+vc3*v3[0] >> vc[1,*]=vc0*v1[1]+vc1*v2[1]+... >> vc[2,*]=vc0*v1[2]+... > >> Ciao, >> Paolo > > well, there's a bit of an index mismatch > in there... > > vc[0,*]=vc0*v1[0]+vc1*v2[0]+vc2*v3[0] > vc[1,*]=vc0*v1[1]+vc1*v2[1]+vc2*v3[1] > vc[2,*]=vc0*v1[2]+vc1*v2[2]+vc2*v3[2] Got it, thanks. I changed it in my test, it turned out to be slower than the original version (about 30% slower). By now, seeing no way to improve the original version with an arbitrary vc, I stick to the version of "adding REBINS" for the case where all indexes from an image are given as input. This is the critical case for my application anyways. When a smaller region is wanted, the original version can be used.
creator: Axel Martínez
date: 2010-09-30T16:38:33-00:00