@kcrisman : That sounds like the right way to do it. Would you by chance be able to point to an example in one of the tutorials that does this in a best-practice way? I tend to confuse myself when converting objects from one class to another, which leads to a lot of overhead, as in my original example. (If you don't have an example at hand, then don't worry, I just had the idea that perhaps you could save me a lot of time)

@kcrisman - thanks for the pointer. Now, most of my code contains much symbolic manipulation before going to numerics, so that's why I tend to start in sage. I will have to do some experimentation ...

Thanks! I guess the take-home-message is to stay in numpy, whereas I was trying to do as much as possible in sage proper, and only shifting to numpy for bottlenecks.

Dear all

I am a sage newbie coming from matlab, R, and Maple, and although I really like many features, I find that my sage code for handling numerical arrays is much more complex and much slower than the corresponding code in R and Matlab.

Here's an example:

Simuating a sample path of standard Brownian motion

I want to simulate a random path of standard Brownian motion. The input is a vector of time points, and I want returned a (random) trajectory, sampled at these time points.

To do this I need: 1) Fast random number generation 2) Fast array operations: Diff, multiplication, and cumsum. I found that in sage I need to go through NumPy, so I wrote:

sage: import numpy
sage: # General routine to sample a Brownian motion on a time grid 
sage: def rB(tvec):
sage:    dt = numpy.diff(numpy.array(tvec))
sage:    sdt = sqrt(dt)
sage:    RNG = RealDistribution('gaussian',1)
sage:    dB = sdt * numpy.array([RNG.get_random_element() for j in range(dt.__len__())],dtype=float)
sage:    B1 = numpy.array([normalvariate(0,sqrt(tvec[0]))],dtype=float)
sage:    B = numpy.cumsum(numpy.concatenate((B1,dB)))
sage:    return(B)

I then time it, testing with 1 million time points:

sage: tvec = srange(0,1e3,1e-3)
sage: time B = rB(tvec)

I get a time of 3.5 seconds on my machine, a standard Laptop running Ubuntu, running sage from the shell. Using prun, I see that half of the time, roughly, is spent generating random numbers and the other half mainly constructing numerical arrays.

For comparison, the R code would be (here I assume that tvec[1]=0)

> tvec <- seq(0,1000,1e-3)
> system.time(B <- c(0,cumsum(rnorm(length(tvec)-1,sd=sqrt(diff(tvec))))))

which executes in 0.24 seconds. The matlab code would be

>> tvec = linspace(0,1e3,1e6+1);
>> tic;B = [0,cumsum(randn(1,length(tvec)-1) .* sqrt(diff(tvec)))];toc

which runs in 0.06 seconds.

You see that my sage code is not only clumsier, it is also between 10 and 100 times slower, so I guess I am doing something wrong.

I would really appreciate if you could point me in the direction of cleaner and faster code.

Cheers,

Uffe