MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Suggestions for improving speed

  • To: mathgroup at smc.vnet.net
  • Subject: [mg113349] Re: Suggestions for improving speed
  • From: Achilleas Lazarides <achilleas.lazarides at gmx.com>
  • Date: Mon, 25 Oct 2010 06:38:51 -0400 (EDT)

Hi,
This seems to be around twice as fast as your version:
===========================================================================
==
ilp == Compile[
   {{j, _Integer}, {jj, _Integer}, {n, _Integer}, {nn, _Integer},
    time, {knots, _Real, 1}, {omyBasis, _Real, 2}},
   Module[
    {jIdx == j + jj,
     myBasis == omyBasis},
    If[nn ==== 1,
     If[knots[[jIdx]] <== time < knots[[jIdx + 1]],
       myBasis[[jj, nn]] == 1;];,
     If[(jj - 1) <== n - (nn - 1),
       myBasis[[jj,
           nn]] == (time - knots[[jIdx]])/(knots[[jIdx + nn - 1]] -
              knots[[jIdx]])*
           myBasis[[jj, (nn - 1)]] + (knots[[jIdx + nn]] -
              time)/(knots[[jIdx + nn]] - knots[[jIdx + 1]])*
           myBasis[[(jj + 1), (nn - 1)]];,
       ];
     ];
    Return[myBasis]
    ]
   ];

NthBasisRecursivenew[{n_, knots_}, j_, time_] :== Module[
  {jj, nn, myBasis},
  myBasis == ConstantArray[0, {n + 1, n + 1}];
  Last@Table[
    Last@Table[
      myBasis == ilp[j, jj, n, nn, time, knots, myBasis],
      {jj, 1, n + 1}
      ],
    {nn, 1, n + 1}
    ];
  myBasis[[1, n + 1]]
  ]
===========================================================================
====

I've just compiled the body of the innermost loop.

Note that changing from your For structures to eg Tables won't help much because most of the time is spent inside the innermost loop anyway (I tried, the change is around 10% faster).

Note also that I don't actually understand what this is supposed to do, nor do I understand what it actually does. So I just literally took what you wrote and moved things around so as to compile. There could be easy ways to make it faster (I don't mean algorithmic, but rather tricks that may or may not apply depending on what is going on).

And finally, all I did was check that it gives a similar-looking plot.

On Oct23, 2010, at 1:04 PM, Joe Hays wrote:

> Hello all,
>
> Today I've been studying B-Splines. I know that Mathematica has a built-in
> function to build the basis of a B-Spline but I wanted to better understand
> it so I coded up my own following the "Cox-de Boor recursion formula" found
> at: http://en.wikipedia.org/wiki/B-spline#Cubic_B-Spline
>
> I tested it against Mathematica's version to make sure things were working
> OK. After I got the right answer I timed the two implementations.
> Mathematica's took ~ 0.12 seconds. My version took ~15.53 seconds. I
> expected my implementation to be slower but not 125X slower!
>
> I realize that I still approach things from the procedural point of view.
> So, I thought I'd ask the Mathematica gurus if they see obvious coding
> choices that cause the drastic slow down...
>
> My tests were:
>
> Given,
>
> m==15;
> n==3;
> knots==Range[0,m-1]
> j==3;
> time==3.5;
>
> Mathematica's version:
> <snip>
> tmpT==AbsoluteTime[];
> Plot[Table[BSplineBasis[{n,knots},j,time],{j,0,m-n-2}],{time,0,m-1},
> PlotRange->All]
> AbsoluteTime[] - tmpT
> </snip>
>
> this took ~0.12 seconds
>
> my version,
> <snip>
> tmpT==AbsoluteTime[];
> Plot[Table[NthBasisRecursive[{n,knots},j,time],{j,0,m-n-2}],{time,0,m-1},
> PlotRange->All]
> AbsoluteTime[] - tmpT
> </snip>
>
> this took ~15.53 seconds
>
> Now, for my code defining "NthBasisRecursive"
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> NthBasisRecursive[{n_, knots_}, j_, time_] :==
> Module[{jj, nn, myBasis, jIdx},
>  myBasis == ConstantArray[0, {n + 1, n + 1}];
>
>  For[nn == 1, nn <== n + 1, nn++,
>   For[jj == 1, jj <== n + 1, jj++,
>
>    jIdx == j + jj; (* Need to offset the knot indexing by the jth basis I'm
> solving for *)
>
>    If[nn ==== 1,
>     If[knots[[jIdx]] <== time < knots[[jIdx + 1]],
>       myBasis[[jj, nn]] == 1;
>       ];
>     ,
>
>     If[(jj - 1) <== n - (nn - 1),
>        myBasis[[jj, nn]] == (time - knots[[jIdx]])/(knots[[jIdx + nn - 1]] -
> knots[[jIdx]])*myBasis[[jj, (nn - 1)]]
>     +  (knots[[jIdx + nn]] -  time)/(knots[[jIdx + nn]] - knots[[jIdx +
> 1]])*myBasis[[(jj + 1), (nn - 1)]];
>      ,
>      ];(* end of inner IF *)
>    ];(* end of outer IF *)
>
>    ];(* end of JJ loop *)
>   ];(* end of nn loop *)
>  myBasis[[1, n + 1]]
>  ]
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> Any insight to help me learn how to make things faster would be greatly
> appreciated!
>
> Joe
>


  • Prev by Date: Re: How to "soft-code" a Block?
  • Next by Date: Re: Can one define own indexing function/mapping to be
  • Previous by thread: Re: Suggestions for improving speed
  • Next by thread: Strange Behaviour of ReplaceAll with Subscript (A bug or what?)