Alec, I was not trying to trash anyone's code or cause offence. I am trying to figure out what was meant by the term "compiler". I see that Fred Lunnon also posted his code here on Mapleprimes, and specifically asked about suggestions for it.

Looking closer, it does not look appropriate to use Compiler:-Compile on the initbagperm procedure. That does indeed set up the arrays (which be switched to Arrays of a suitable datatype) which then get passed into procedure nextbagperm which does the hard work. So nextbagperm would be the natural Compiler target.

Initial experimentation indicates some speed-up buy using Compiler:-Compile. At least a factor of 5 at a smaller size. I've posted that as a reply to the original code as posted here on Mapleprimes.

acer

One would want to double-check the actual results, to make sure that each version of the code below produced the same thing.

Here is a run of three versions: the original nextbagperm, a "typed" nextbagperm, and a Maple-Compiled nextbagperm. The second is about seven times slower than the first. The third is about five times faster than the first.

    |\^/|     Maple 12 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2008
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> initbagperm := proc (rank)
> global  R_, P_, Q_, D_;
> local m, j, r, r_1, r_m;
> m := nops(rank); r := sort(rank);
> if m <= 0 then print(`bagperm: length <= 0`)
> else r_1, r_m := r[1], r[m];
> R_ := Array(0..m+3, [r_1-1, seq(r[j], j = 1..m), r_1-1, r_m+1, r_m+2],datatype=integer[4]);
> P_ := Array(0..m+3, [0, seq(j, j = 1..m), 0, m+1, m+2],datatype=integer[4]); # permed identities
> Q_ := Array(0..m+3, [0, seq(j, j = 1..m), m+2, m+3, 0],datatype=integer[4]); # inverse locations
> D_ := Array(0..m+3, [0, seq(+1, j = 1..m), 0, +1, +1],datatype=integer[4]); # permed drifts
> fi; NULL end:
>
> nextbagperm := proc (R_, P_, Q_, D_)
> local i, j, k, l, d, e, p;
# locate earliest unblocked element at j, starting at blocked element 0
> j, i, d := 0, 0, 0;
> while R_[j] >= R_[i] do
> D_[j] := -d; # blocked at j; reverse drift d pre-emptively
> j := Q_[P_[j]+1]; d := D_[j]; i := j+d; # next element at j, neighbour at i
> if R_[j-1] <> R_[j] then l := j # save left end of run in l
> elif d < 0 then i := l-1 fi od; # restore left end at head of run
# shift run of equal rank from i-d,i-2d,...,l to i,i-d,...,l+d
> k := i; if d < 0 then l := j fi;
> e := D_[i]; p := P_[i]; # save neighbour drift e and identifier p
> while k <> l do
> P_[k] := P_[k-d]; Q_[P_[k]] := k;
> D_[k] := -1; k := k-d od; # reset drifts of run tail elements
> R_[l], R_[i] := R_[i], R_[l]; # transpose user ranks
> D_[l], D_[i] := e, d; # restore drifts of head and neighbour
> P_[l], Q_[p] := p, l; # wrap neighbour around to other end
> j end:
>
> nextbagperm_typed := proc (R_::Array(0..18,datatype=integer[4]),
>                      P_::Array(0..18,datatype=integer[4]),
>                      Q_::Array(0..18,datatype=integer[4]),
>                      D_::Array(0..18,datatype=integer[4]))
> local i :: integer, j :: integer, k :: integer, l :: integer,
>       d :: integer, e :: integer, p :: integer;
# locate earliest unblocked element at j, starting at blocked element 0
> j, i, d := 0, 0, 0;
> while R_[j] >= R_[i] do
> D_[j] := -d; # blocked at j; reverse drift d pre-emptively
> j := Q_[P_[j]+1]; d := D_[j]; i := j+d; # next element at j, neighbour at i
> if R_[j-1] <> R_[j] then l := j # save left end of run in l
> elif d < 0 then i := l-1 fi od; # restore left end at head of run
# shift run of equal rank from i-d,i-2d,...,l to i,i-d,...,l+d
> k := i; if d < 0 then l := j fi;
> e := D_[i]; p := P_[i]; # save neighbour drift e and identifier p
> while k <> l do
> P_[k] := P_[k-d]; Q_[P_[k]] := k;
> D_[k] := -1; k := k-d od; # reset drifts of run tail elements
> R_[l], R_[i] := R_[i], R_[l]; # transpose user ranks
> D_[l], D_[i] := e, d; # restore drifts of head and neighbour
> P_[l], Q_[p] := p, l; # wrap neighbour around to other end
> j end:
>
> compiled_nextbagperm_typed := Compiler:-Compile(nextbagperm_typed):
>
> parts := [5, 5, 5]; # chosen partition of set
                              parts := [5, 5, 5]
 
> n := nops(parts): m := add(parts[i], i = 1..n): # length and weight
> symb := [seq(seq(i, j = 1..parts[i]), i = 1..n)]:
> cardi := combinat[multinomial](m, op(parts)): # number of bag-perms
>
> kernelopts(printbytes=false);
                                     true
 
>
> st:=time():
> initbagperm(symb, R_,P_,Q_,D_):
> n := 1: print(n, [seq(R_[i], i = 1..m)]):
               1, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
 
> while nextbagperm(R_,P_,Q_,D_) <= m and n < cardi + 2 do
>   n := n+1;
>   #print(n, [seq(R_[i], i = 1..m)]);
> od:
> print(n, cardi); # equal
                                756756, 756756
 
> time()-st;
                                    24.098
 
>
> st:=time():
> initbagperm(symb, R_,P_,Q_,D_):
> n := 1: print(n, [seq(R_[i], i = 1..m)]);
               1, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
 
> while nextbagperm_typed(R_,P_,Q_,D_) <= m and n < cardi + 2 do
>   n := n+1;
>   #print(n, [seq(R_[i], i = 1..m)]);
> od:
> print(n, cardi); # equal
                                756756, 756756
 
> time()-st;
                                    173.600
 
>
> st:=time():
> initbagperm(symb, R_,P_,Q_,D_):
> n := 1: print(n, [seq(R_[i], i = 1..m)]);
               1, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
 
> while compiled_nextbagperm_typed(R_,P_,Q_,D_) <= m and n < cardi + 2 do
>   n := n+1;
>   #print(n, [seq(R_[i], i = 1..m)]);
> od:
> print(n, cardi); # equal
                                756756, 756756
 
> time()-st;
                                     5.239

As one can see, I had a little difficulty with the sizes in the type-specification of the Array arguments to the nextbagperm procedure. If I did not hard code the size then Compile complained that the lower index started at 0 instead of 1. It may be that the code would need rewriting to be "1-based" instead of "0-based" for Array referencing, so that the Array size need not be supplied in the parameter type specification. I'd welcome any news that there is another way.

acer

One would want to double-check the actual results, to make sure that each version of the code below produced the same thing.

Here is a run of three versions: the original nextbagperm, a "typed" nextbagperm, and a Maple-Compiled nextbagperm. The second is about seven times slower than the first. The third is about five times faster than the first.

    |\^/|     Maple 12 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2008
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> initbagperm := proc (rank)
> global  R_, P_, Q_, D_;
> local m, j, r, r_1, r_m;
> m := nops(rank); r := sort(rank);
> if m <= 0 then print(`bagperm: length <= 0`)
> else r_1, r_m := r[1], r[m];
> R_ := Array(0..m+3, [r_1-1, seq(r[j], j = 1..m), r_1-1, r_m+1, r_m+2],datatype=integer[4]);
> P_ := Array(0..m+3, [0, seq(j, j = 1..m), 0, m+1, m+2],datatype=integer[4]); # permed identities
> Q_ := Array(0..m+3, [0, seq(j, j = 1..m), m+2, m+3, 0],datatype=integer[4]); # inverse locations
> D_ := Array(0..m+3, [0, seq(+1, j = 1..m), 0, +1, +1],datatype=integer[4]); # permed drifts
> fi; NULL end:
>
> nextbagperm := proc (R_, P_, Q_, D_)
> local i, j, k, l, d, e, p;
# locate earliest unblocked element at j, starting at blocked element 0
> j, i, d := 0, 0, 0;
> while R_[j] >= R_[i] do
> D_[j] := -d; # blocked at j; reverse drift d pre-emptively
> j := Q_[P_[j]+1]; d := D_[j]; i := j+d; # next element at j, neighbour at i
> if R_[j-1] <> R_[j] then l := j # save left end of run in l
> elif d < 0 then i := l-1 fi od; # restore left end at head of run
# shift run of equal rank from i-d,i-2d,...,l to i,i-d,...,l+d
> k := i; if d < 0 then l := j fi;
> e := D_[i]; p := P_[i]; # save neighbour drift e and identifier p
> while k <> l do
> P_[k] := P_[k-d]; Q_[P_[k]] := k;
> D_[k] := -1; k := k-d od; # reset drifts of run tail elements
> R_[l], R_[i] := R_[i], R_[l]; # transpose user ranks
> D_[l], D_[i] := e, d; # restore drifts of head and neighbour
> P_[l], Q_[p] := p, l; # wrap neighbour around to other end
> j end:
>
> nextbagperm_typed := proc (R_::Array(0..18,datatype=integer[4]),
>                      P_::Array(0..18,datatype=integer[4]),
>                      Q_::Array(0..18,datatype=integer[4]),
>                      D_::Array(0..18,datatype=integer[4]))
> local i :: integer, j :: integer, k :: integer, l :: integer,
>       d :: integer, e :: integer, p :: integer;
# locate earliest unblocked element at j, starting at blocked element 0
> j, i, d := 0, 0, 0;
> while R_[j] >= R_[i] do
> D_[j] := -d; # blocked at j; reverse drift d pre-emptively
> j := Q_[P_[j]+1]; d := D_[j]; i := j+d; # next element at j, neighbour at i
> if R_[j-1] <> R_[j] then l := j # save left end of run in l
> elif d < 0 then i := l-1 fi od; # restore left end at head of run
# shift run of equal rank from i-d,i-2d,...,l to i,i-d,...,l+d
> k := i; if d < 0 then l := j fi;
> e := D_[i]; p := P_[i]; # save neighbour drift e and identifier p
> while k <> l do
> P_[k] := P_[k-d]; Q_[P_[k]] := k;
> D_[k] := -1; k := k-d od; # reset drifts of run tail elements
> R_[l], R_[i] := R_[i], R_[l]; # transpose user ranks
> D_[l], D_[i] := e, d; # restore drifts of head and neighbour
> P_[l], Q_[p] := p, l; # wrap neighbour around to other end
> j end:
>
> compiled_nextbagperm_typed := Compiler:-Compile(nextbagperm_typed):
>
> parts := [5, 5, 5]; # chosen partition of set
                              parts := [5, 5, 5]
 
> n := nops(parts): m := add(parts[i], i = 1..n): # length and weight
> symb := [seq(seq(i, j = 1..parts[i]), i = 1..n)]:
> cardi := combinat[multinomial](m, op(parts)): # number of bag-perms
>
> kernelopts(printbytes=false);
                                     true
 
>
> st:=time():
> initbagperm(symb, R_,P_,Q_,D_):
> n := 1: print(n, [seq(R_[i], i = 1..m)]):
               1, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
 
> while nextbagperm(R_,P_,Q_,D_) <= m and n < cardi + 2 do
>   n := n+1;
>   #print(n, [seq(R_[i], i = 1..m)]);
> od:
> print(n, cardi); # equal
                                756756, 756756
 
> time()-st;
                                    24.098
 
>
> st:=time():
> initbagperm(symb, R_,P_,Q_,D_):
> n := 1: print(n, [seq(R_[i], i = 1..m)]);
               1, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
 
> while nextbagperm_typed(R_,P_,Q_,D_) <= m and n < cardi + 2 do
>   n := n+1;
>   #print(n, [seq(R_[i], i = 1..m)]);
> od:
> print(n, cardi); # equal
                                756756, 756756
 
> time()-st;
                                    173.600
 
>
> st:=time():
> initbagperm(symb, R_,P_,Q_,D_):
> n := 1: print(n, [seq(R_[i], i = 1..m)]);
               1, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
 
> while compiled_nextbagperm_typed(R_,P_,Q_,D_) <= m and n < cardi + 2 do
>   n := n+1;
>   #print(n, [seq(R_[i], i = 1..m)]);
> od:
> print(n, cardi); # equal
                                756756, 756756
 
> time()-st;
                                     5.239

As one can see, I had a little difficulty with the sizes in the type-specification of the Array arguments to the nextbagperm procedure. If I did not hard code the size then Compile complained that the lower index started at 0 instead of 1. It may be that the code would need rewriting to be "1-based" instead of "0-based" for Array referencing, so that the Array size need not be supplied in the parameter type specification. I'd welcome any news that there is another way.

acer

This order of this thread is a little mixed up. But this comment above may help.

acer

This order of this thread is a little mixed up. But this comment above may help.

acer

What is the "compiler" that gets mentioned on that page? Maple uses an interpreter, and doesn't (normally) compile code. Since the usual (not evalhf) Maple interpreter doesn't compile (not even just in time), it doesn't seem surprising to me that adding type-checks incurs extra runtime cost. It seems like a misunderstanding about how Maple works. That's not any disparagement -- Fred Lunnon did not create Maple in the way that Knuth wrote TeX.

Maple did get a Compiler routine, a few years ago. But Compiler:-Compile doesn't know how to handle that code, as it stands. (The Compiler knows about Arrays, but not arrays. And it might be more feasible to have the procedure accept the Arrays as arguments, since the Compiler doesn't know how to handle indeterminate index ranges like 0..m+3.)

With a little work, an edited version of the procedure (possibly running inside a wrapping procedure that sets up the Arrays) might be done. I would expect that to get a speed-up of something like a factor of 10-25 or so, for the sub-parts that are "Compilable". One would have to hit it explicitly with Compiler:-Compile.

acer

If I changed elements like `-p/2` to `- p/2`, with a space after the minus sign, then it appeared better in the exported .eps file. The symbol showed as minus rather than plusminus, and the extra whitespace didn't seem too glaring.

But then I tried something else. There is also a programmatic way to export the plot. before the display call, I added this command,

plotsetup(ps,plotoutput=`foo.eps`,plotoptions=`portrait,noborder`);

That command will make plots get written directly to the file. The plot won't also appear visually. Using the exporter, I found that the original syntax you used of `-p/2` would come out with the correct minus sign in the .eps file.

I used ghostview to preview the .eps file, rather than embedding it into latex.

Is the "90 degree rotation" issue just one of portrait versus landscape, for Postscript?

It's fun, no, that Maple has (at least) three distinct Postscript converters? (Ie, the menu driven one in Standard, the menu driven one in Classic, and that programmatic command available in each interface.)

The command plotsetup(default) will reset the plot device to its original default setting (ie, visual display of plots).

acer

If I changed elements like `-p/2` to `- p/2`, with a space after the minus sign, then it appeared better in the exported .eps file. The symbol showed as minus rather than plusminus, and the extra whitespace didn't seem too glaring.

But then I tried something else. There is also a programmatic way to export the plot. before the display call, I added this command,

plotsetup(ps,plotoutput=`foo.eps`,plotoptions=`portrait,noborder`);

That command will make plots get written directly to the file. The plot won't also appear visually. Using the exporter, I found that the original syntax you used of `-p/2` would come out with the correct minus sign in the .eps file.

I used ghostview to preview the .eps file, rather than embedding it into latex.

Is the "90 degree rotation" issue just one of portrait versus landscape, for Postscript?

It's fun, no, that Maple has (at least) three distinct Postscript converters? (Ie, the menu driven one in Standard, the menu driven one in Classic, and that programmatic command available in each interface.)

The command plotsetup(default) will reset the plot device to its original default setting (ie, visual display of plots).

acer

Users will expect that their old code works, as is. And the syntax of the commands from the two packages is different. So the direct substitution of the package names, in commands, would not work I think.

One can observe that effort has been made to more strongly deprecate linalg, in the Maple help-system in Maple 12. For example, entering,

> ?linalg

> ?stats

takes one to the Student[LinearAlgebra] and Statistics help-pages, repectively, in Maple 12.

It's a little heavy handed in parts. For example issuing ?array takes one to the Array help-page, but I don't see any obvious cross-link from the top of that page to ?array(deprecated). One can find similar problems, when trying to access other older deprecates pages in Maple 12 without knowing in advance that they get accessed by the syntax ?<fcn>(deprecated). In the Standard interface, the help-browser should show a list of relevant matches, but in the Classic and commandline interfaces it's more problematic.

acer

I suspect that he is referring to the Mapleprimes points system.

For example, by clicking on the number beside your "red leaf". And so one can get to the rankings here.

To try to answer the question, about what they are used for, not much except to indicate how many posts have been made from the given member-handle. I believe that some ability to edit own's own posts (or perhaps to post blog entries?) is enabled once one attains a number (10?) of posts.

acer

I suspect that he is referring to the Mapleprimes points system.

For example, by clicking on the number beside your "red leaf". And so one can get to the rankings here.

To try to answer the question, about what they are used for, not much except to indicate how many posts have been made from the given member-handle. I believe that some ability to edit own's own posts (or perhaps to post blog entries?) is enabled once one attains a number (10?) of posts.

acer

One of the things that makes textbooks on Maple tricky is that the system keeps evolving. The linalg package gets replaced with LinearAlgebra, the stats package with Statistics, the networks package with GraphTheory, Maplets with embedded components, etc. And some parts get enhanced, for example the newer DAE solvers. So, given a mathematical or programming task, it's quite conceivable that the nicer ways to do it in Maple may change with the advent of a new release. And Maple has been getting new releases most every year, for some time.

So, a text written in 2002 (like the one you cited) may be partially out of date. I find it hard to tell whether the more general Maple books or the ones focussed on some discipline fall behind current maple functionality more quickly. The title of the book that you cited mentions probability and stochastic DEs. The newer Statistics package postdates the book, I believe (as does the Financial Modelling Toolbox, which does Wiener process stuff, etc).

So, depending on the hints you get (here, say) about changes in the product and the books age, you may wish to be extra careful in judging a particular book.

acer

One of the things that makes textbooks on Maple tricky is that the system keeps evolving. The linalg package gets replaced with LinearAlgebra, the stats package with Statistics, the networks package with GraphTheory, Maplets with embedded components, etc. And some parts get enhanced, for example the newer DAE solvers. So, given a mathematical or programming task, it's quite conceivable that the nicer ways to do it in Maple may change with the advent of a new release. And Maple has been getting new releases most every year, for some time.

So, a text written in 2002 (like the one you cited) may be partially out of date. I find it hard to tell whether the more general Maple books or the ones focussed on some discipline fall behind current maple functionality more quickly. The title of the book that you cited mentions probability and stochastic DEs. The newer Statistics package postdates the book, I believe (as does the Financial Modelling Toolbox, which does Wiener process stuff, etc).

So, depending on the hints you get (here, say) about changes in the product and the books age, you may wish to be extra careful in judging a particular book.

acer

I would like to thank all who posted to this thread. The points raised are very interesting. Hopefully, having several major contributors to this site offer their suggestions/needs will induce those ideas to be weighed and judged carefully.

ps. I'm still editing my list...

acer

LinearSolve uses a general (full rectangular dense) solver on tridiagonal linear floating-point systems. For real data, that means LAPACK dgetrf and dgetrs (NAG f07adf and f07aef). But LinearSolve could be taught to instead use dgttrf and dgttrs which are specialized.

> infolevel[LinearAlgebra]:=1:
> with(LinearAlgebra):

> f := n->Matrix(n,n,shape=band[1,1],scan=band[1,1],
>    [[1$(n-1)],[seq(i,i=1..n)],[1$(n-1)]],
>    storage=rectangular,datatype=float[8]):

> rtable_indfns(f(5)); # it's tridiagonal in structure...
                                  band[1, 1]
 
> rtable_options(f(5),storage); # ...and full storage, as dgttrf expects
                                  rectangular

> LinearSolve(evalf(f(5)),evalf(Vector(5,i->i)));
LinearSolve:   "using method"   LU
LinearSolve:   "calling external function"
LinearSolve:   "NAG"   hw_f07adf
LinearSolve:   "NAG"   hw_f07aef
                            [0.696969696969697017]
                            [                    ]
                            [0.303030303030302983]
                            [                    ]
                            [0.696969696969697017]
                            [                    ]
                            [0.606060606060606077]
                            [                    ]
                            [0.878787878787878896]

Moreover, if the (real) Matrix is symmetric positive-definite then dpttrf and dpttrs could be used.

> type(f(5),'Matrix'(symmetric));
                                     true

> IsDefinite(evalf(f(5)),query=positive_definite);
                                     true

acer