The resulting matrix is symmetric in the back diagonal.
This (rough) solution generates the upper triangle, much faster and with a
larger footprint.
cAO=: 4 : 0
z=.(|:x#,:t),:(t=.i.#y)+/|.i.x
v=. */(((<:x)#0),~ t){~ z
w=. s%~;(>:i.x)<@{."_1 &.|:(s=.>:(#y)-x)+/\ v
(2#x)$(w,((-:@*<:)x)#0)(;</.r)},r=.i.2#x
)
ts'K=:276 ([: ((+ (*/ |.))~/ % <:@#) ]\ , 0:) i.600'
0.51149222 2108992
K-: |.@|:@|.K NB. symmetric in back diagonal
1
(276 cAO i.600) -: ((]>[:+/~i.)@#*]) K NB. upper triangle
of K
1
ts'276 cAO i.600'
0.024038428 10516928
0.51149222 2108992%0.024038428 10516928
21.278106 0.20053308
ts'K=:276 ([: ((+ (*/ |.))~/ % <:@#) ]\ , 0:) i.4000'
6.1219979 8408576
ts'276 cAO i.4000'
0.11301888 58776768
(276 cAO i.4000)-:((]>[:(+)/~i.)@#*])K
1
Overall efficiency improvement:
*/ 6.1219979 8408576 % 0.11301888 58776768
7.7492368
R.E. Boss
> -----Oorspronkelijk bericht-----
> Van: [EMAIL PROTECTED] [mailto:programming-
> [EMAIL PROTECTED] Namens Oleg Kobchenko
> Verzonden: donderdag 8 mei 2008 9:25
> Aan: Programming forum
> Onderwerp: Re: [Jprogramming] Outer products problem
>
> Here's one without boxing.
> The idea of using zero element to defer
> expansion to run is great.
>
> ts '276 calcAvgOP i.600'
> 0.722269 1.35284e8
>
> ts '276 calcAvgOP2 i.600'
> 0.534168 1.71264e6
>
> ts '325 %~ > ((*/|.)@[ + ])&.> / (276 <\ i.600),<0'
> 0.528031 2.25075e6
>
> ts '276 ([: ((+ (*/ |.))~/ % <:@#) ]\ , 0:) i.600' NB. <--
> 0.529847 2.10893e6
>
>
> --- On Wed, 5/7/08, Roger Hui <[EMAIL PROTECTED]> wrote:
>
> > ts '276 calcAvgOP i.600'
> > 1.1417 1.35284e8
> > ts '325 %~ > ((*/|.)@[ + ])&.> / (276
> > <\ i.600),<0'
> > 0.734065 2.25075e6
> >
> >
> >
> > ----- Original Message -----
> > From: Ronan Reilly <[EMAIL PROTECTED]>
> > Date: Wednesday, May 7, 2008 12:02
> > Subject: [Jprogramming] Outer products problem
> > To: Programming forum <[email protected]>
> >
> > > I'm doing some analyses on EEG data and have run
> > into a limit
> > > error for
> > > vectors of > 4000 elements. I'm wondering if
> > there's any
> > > way around it.
> > >
> > > The code involves moving a fixed width window across a
> > long
> > > vector of values
> > > one element at a time, computing the outer product of
> > the
> > > (reversed) window
> > > of values with itself, and averaging over all the
> > resulting
> > > matrices. I run
> > > into the limit error for vectors greater than 4000
> > elements with
> > > a typical
> > > window size of 276 - yes, I was naive not to expect
> > this :-)
> > >
> > > Here's the code that hits the limit:
> > >
> > > NB. x = window size ; y = vector of EEG values
> > > calcAvgOP =: [: (+/ % #) (*"0 1 |.)\
> > >
> > > $M =: 276 calcAvgOP i. 4000
> > > |limit error: calcAvgOP
> > > | M=:276 calcAvgOP i.4000
> > >
> > > The problem is having to hold onto the millions of
> > intervening
> > > matricesbefore averaging. I'd very much like to
> > find a
> > > tacit solution to this
> > > problem and I was wondering if there's a
> > formulation that can
> > > circumvent it.
> > >
> > > On the other hand, the rather inelegant explicit
> > definition
> > > below works
> > > fine.
> > >
> > > NB. x = window size ; y = vector of EEG values
> > > calcAvgOP2 =: dyad define
> > > i =. 0
> > > lim =. >: x -~ $y
> > > acc =. 0 $~ x,x
> > > while. i < lim do.
> > > acc =. acc + (*"0 1
> > > |.) x {. y
> > > y =. }. y
> > > i =. >: i
> > > end.
> > > acc % i
> > > )
> > >
> > > $M =: 276 calcAvgOP2 i. 4000
> > > 276 276
> > ----------------------------------------------------------------------
> > For information about J forums see
> > http://www.jsoftware.com/forums.htm
>
>
>
> __________________________________________________________________________
> __________
> Be a better friend, newshound, and
> know-it-all with Yahoo! Mobile. Try it now.
> http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
