Few questions (convolutions, 'mapping', linspace)

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Few questions (convolutions, 'mapping', linspace)

Przemek Klosowski-3
Hi,

I have a couple of questions that came in my recent using of Octave; I wonder
if anyone could help:

 - I want to convolute the theoretical function S(Q,w) with the experimental
   resolution r(w), to compare the resulting Sr(Q,w) with experimental data.
   The formula is:
                  w0
                  /
      Sr(Q,w) =   |  S(Q,w-v) r(v) dv
                  /
                 -w0

   I was hoping to use the quad("fun",-w0,w0), where fun would be defined
   as 'function y=fun(x); y=S(Q,w-x)*r(x); endfunction'. This didn't quite
   work, as I need to access Q,w from within that function. I can't pass
   them as parameters of fun() because quad() expects only one-parameter fun;
   I couldn't figure the way to pass them as globals either.

 - I often find myself defining a function (y=fun(x)), and hoping to be able
   to 'map' (in the LISP sense) the vectors/matrices with this function, for
   instance to plot it. The naive approach, y=f(x) for a vector/matrix 'x',
   does not result in

         y = [f(x(1)),  f(x(2)),  f(x(3)) ... f(x(N))],

   which is fine because the arithmetic operators in the function definition
   have meaning for matrices as well. For vector 'x' I can always do

         for i=1:length(x); y(i)=f(x(i)); endfor

   but I wonder if there is a neater, general way. I suppose I could be defining
   my functions to do the explicitly recognize vector arguments, but
   that still leaves the question whether I have to write explicit 'for' loops
   or is there a better way.
   
 - I noticed that the linspace operator is fairly slow: linspace(-1,1, 10000) takes
   25 seconds on a 150MHz r4400, while the a=(0:.0001:1); is very quick, and the
   FFT of the resulting matrix is only few seconds. What is the advantage of linspace
   over the range operator, and why is it so slow?



        Greetings
                        przemek

Reply | Threaded
Open this post in threaded view
|

Re: Few questions (convolutions, 'mapping', linspace)

John Eaton-4
[hidden email] (Przemek Klosowski) wrote:

:  - I noticed that the linspace operator is fairly slow:
:    linspace(-1,1, 10000) takes 25 seconds on a 150MHz r4400, while
:    the a=(0:.0001:1); is very quick, and the FFT of the resulting
:    matrix is only few seconds. What is the advantage of linspace  
:    over the range operator, and why is it so slow?

The range operator is built-in, and in Octave it creates a special
range object that only contains the base, limit, and increment, so it
requires the same amount of work and storage to say

  1:3:5

or

  1:0.1:1e6

The linspace function is currently implemented as an external
function, and creates a matrix with a for loop:

    delta = (x2 - x1) / (npoints - 1);
    retval = zeros (1, npoints);
    for i = 0:npoints-1
      retval (i+1) = x1 + i * delta;
    endfor

which is terribly inefficient.  Changing it to

    delta = (x2 - x1) / (npoints - 1);
    retval = x1:delta:x2;

should not change the result and will speed things up considerably.

With this change, about the only advantage of the linspace function is
that you don't have to compute the increment yourself.

BTW, the loop in logspace() should probably also be replaced with

    retval = 10 .^ retval;

Thanks for pointing this out.

jwe

Reply | Threaded
Open this post in threaded view
|

Re: Few questions (convolutions, 'mapping', linspace)

John Eaton-4
In reply to this post by Przemek Klosowski-3
[hidden email] (Przemek Klosowski) wrote:

:  - I often find myself defining a function (y=fun(x)), and hoping to
:    be able to 'map' (in the LISP sense) the vectors/matrices with
:    this function, for instance to plot it.

Here's one way to do it, but it won't be very fast.

  function y = map (f, x)

  # Given a function y = f (x) that expects a scalar argument and
  # produces a scalar result, call it for each element of the matrix x
  # and return the corresponding matrix of values.

    if (nargin != 2)
      error ("usage: map (f, x)");
    endif

    if (! isstr (f))
      error ("map: expecting string as first argument");
    endif

    [nr, nc] = size (x);
    y = zeros (nr, nc);
    for j = 1:nc
      for i = 1:nr
        y (i, j) = feval (f, x (i, j));
      endfor
    endfor

  endfunction

It might be desirable to add this as a built-in function to speed it
up, but even then it won't be all that fast because it will require
multiple calls to a user-supplied function, which is fairly slow.

jwe

Reply | Threaded
Open this post in threaded view
|

Re: Few questions (convolutions, 'mapping', linspace)

John Eaton-4
In reply to this post by Przemek Klosowski-3
[hidden email] (Przemek Klosowski) wrote:

:  - I want to convolute the theoretical function S(Q,w) with the experimental
:    resolution r(w), to compare the resulting Sr(Q,w) with experimental data.
:    The formula is:
:                   w0
:                   /
:       Sr(Q,w) =   |  S(Q,w-v) r(v) dv
:                   /
:                  -w0
:
:    I was hoping to use the quad("fun",-w0,w0), where fun would be defined
:    as 'function y=fun(x); y=S(Q,w-x)*r(x); endfunction'. This didn't quite
:    work, as I need to access Q,w from within that function. I can't pass
:    them as parameters of fun() because quad() expects only one-parameter fun;
:    I couldn't figure the way to pass them as globals either.

It should work to write

  function y = S (Q, w, v)
   y = ...
  endfunction

  function y = r (v)
   y = ...
  endfunction

  global Q;
  global w;

  Q = ...
  w = ...

  function y = fun (x)
    global Q;
    global w;
    y = S (Q, w-x) * r (x);
  endfunction

  w0 = ...

  quad ("fun", -w0, w0)


If this fails, please send a complete bug report to bug-octave.

Thanks,

jwe