Re: Convolution of Stick Spectra

[mole6e23(at)hotmail.com (Todd Clements)]

craigmnet@cow.physics.wisc.edu wrote:

>After looking at your problem it looks like the widths of the lines
>are as wide as your energy range.  I guess this would be appropriate
>[snip..]
>so many exponentiations.  You could have gotten a pretty big savings
>if the lines were narrow, and you could restrict the computation to a
>narrow region around the line center (say +/- 10 sigma with WHERE).  I
>recommend that anyway to get rid of underflow errors.
>[snip..]
>If for example your sigma term were
>(((.12*sqrt(stick[0,i]/1000))/1.6651)*10) [ note the last factor is
>smaller ] then things start to look interesting.

That's because there's always a danger in taking code x and modifying it
into simpler example y! As you suggested, I actually DO have a sigma term
with a *10 instead of a *1000.

I did what you suggested with the +/- 10*sigma with where (code below),
and it drastically improved running time. For the 5000 element array, the
running time went from 62.3 seconds to 4.9 seconds! I plot the error
associated with the method, and I was actually able to go down to +/-
4*sigma before there was any noticable error.

Thanks!
Todd

--

pro test

 convoluted = fltarr( 2, 2048 )
 convoluted[0,*] = findgen( 2048 ) / 2047. * 4.
 convoluted2 = convoluted
 
 ;; Let's fake a stick spectrum, we usually have at least this many elements
 stick = abs(randomn(10, 2, 5000 ))
 stick[0,*] = stick[0,*] * 4.
 stick[1,*] = stick[1,*] * 1000.

 time=systime(1) 
 for i=0L, n_elements( stick ) / 2 -1 do $
      convoluted[1,*] = stick[1,i]*exp(-((convoluted[0,*] - stick[0,i])^2)/ $
                  (((.12*sqrt(stick[0,i]/1000))/1.6651)*10)^2) $
                  + convoluted[1,*]

 time2 = systime(1)
 for i=0L, n_elements( stick ) / 2 -1 do begin
     ;; Find range
     sigma = (((.12*sqrt(stick[0,i]/1000))/1.6651)*10)
     range = where( convoluted2[0,*] ge (stick[0,i]-4*sigma) and $
                    convoluted2[0,*] le (stick[0,i]+4*sigma) )

     if( range[0] ne -1 ) then $
       convoluted2[1,range] = stick[1,i]*exp(-((convoluted2[0,range] -  $
                             stick[0,i])^2)/sigma^2) $
                             + convoluted2[1,range]
 endfor

 print, 'time old: ', time2-time
 print, 'time new: ', systime(1)-time2
 
 new = convoluted
 new[1,*] = new[1,*] - convoluted2[1,*]

 plot,new[0,*], new[1,*],yrange=[-10,10]

end ;; test