Tools of the Maple Masters

Click "Add child page" to create a new item.

building up sets and removing redundant elements

"I've seen this element before..."

Often we are faced with the problem of building up sets incrementally, by removing pieces one at a time from a larger whole. The bottlenecks in this case are usually:
1) adding a small set X to a large set S (copies S and X, making this ~O(|S|+|X|))
2) removing elements of the large set S from the small set X (binary search: |X|*log(|S|))

A classic example of this is a breadth-first-search. We start at one vertex of a graph and in each iteration we add the set of new neighbors X to the set of vertices S that have already been found. We can make this more useful by making the program return the sets of new neighbors found in each iteration, that is, the sets of vertices that are distance 1, 2, 3, etc. from the initial vertex.

Here is a simple example of a graph:

1--2--3
|     |   
4--5--6
   |  |
7--8--9

which I will store as a table mapping vertices (e.g.: 1) to a set of neighbors (e.g.: {2,4}).

G := table([1={2,4}, 2={1,3}, 3={2,6}, 4={1,5}, 5={4,6,8}, 6={3,5,9}, 7={8}, 8={5,7,9}, 9={6,8}]);

Starting from vertex 1, the BFS should return [{1}, {2,4}, {3,5}, {6,8}, {7,9}], corresponding to vertices that are distance [0,1,2,3,4] away.

How would you program this in Maple ? I like to use the following trick:

zap := proc(v) procname(v) := NULL; v end proc:
forget(zap):

This procedure returns any object that it has not seen before, but assigns object=NULL to its own remember table so that if the object is encountered again nothing is returned. If S is a set, the result of map(zap, S) is the set of objects in S that are new. If you use this in your code, you should always take care to include the second line, forget(zap). Otherwise you will find yourself tracking down mysterious bugs, the cause of which are a topic in itself.

As an aside, what useful operation will the following do ?

zap := proc(v) procname(v) := NULL; v end proc:
forget(zap):
L := [1,2,1,1,3,2,4,3,1,2];
map(zap, L);

And for the experts in the house, how might you use the following variant ?

zap := proc(v) procname(v) := false; true end proc:
forget(zap);

Returning to breadth-first-search, here is some simple yet efficient code that takes a graph G as encoded above and a vertex v to start from. We store the current new neighbors in a set X and write them to a table S.

BFS := proc(G, v)
  local S, X, i, zap, d;
   zap := proc(v) procname(v) := NULL; v end proc;
   forget(zap);
   S := table();
   X := map(zap, {v});
   for i from 0 while nops(X) > 0 do
     S[i] := X;
     X := {seq(op(map(zap, G[i])), i=X)};
   end do;
   d := i-1;
   [seq(S[i], i=0..d)];
end proc:

You can play with this code to make it more efficient. For example, you can use another seq in place of op(map(zap, G[i])) and use lists everywhere instead of sets to completely avoid sorting (all elements are unique). These types of optimizations are important for large problems. Assuming you know, you should always ask yourself "what will the kernel do ?" That's the surest way to fast code.

Cartesian Products of Lists

The Cartesian product of a sequence of m lists L₁, ..., L_n, with list i having n_i elements, is given by the sequence of Π_i n_i lists [L₁₁,...,L_m1], ..., [L_1n₁,...,L_{mn_m}], where L_jk is the k-th element of the j-th list.

If the right-most element in the lists varies the fastest, the sequence is in row-major order, if the left-most elements varies the fastest, the sequence is in column-major order. With one exception (noted) all the routines given below use row-major order.

Iterators

combinat[cartprod]

Because a Cartesian product can be large, Maple provides a procedure, combinat[cartprod], that generates an iterator for the Cartesian product. Each call to this iterator returns the next value in the sequence. Here is a procedure that uses cartprod to iterate through each list in the Cartesian product of its arguments, passing each generated list to the function f.

Note: This procedure, as well as most of the following, uses the seq parameter modifier, a mechanism introduced in Maple 11. The parameter L is assigned an expression sequence consisting of the following lists. Here it is equivalent to using the older args[2..-1], provided the remaining arguments are lists.

use_combinat := proc(f, L::seq(list))
local T;
uses combinat;
    T := cartprod([L]);
    while not T['finished'] do
        f(T['nextvalue']());
    end do;
    NULL;
end proc:

By assigning f to print its argument, we can demonstrate that this works as intended:

f := curry(printf,"  %a"):
lsts := [$1..3],[a,b]:
use_combinat(f,lsts): printf("\n"):
  [1, a]  [1, b]  [2, a]  [2, b]  [3, a]  [3, b]

Using combinat[cartprod] is a straight-forward and effective way to step through the entire set without having to allocate memory for the set. The iterator itself, however, requires some time to execute and if the loop itself is otherwise fast—say it calls a compiled procedure—then the time spent evaluating the iterator can be significant.

rtable_scanblock

The Maple built-in procedure rtable_scanblock applies a function to each index of an rtable. If each dimension of the rtable is a range that corresponds to the number of elements in each of the lists, then there is a simple mapping from the indices of the rtable to the sublists in the Cartesian product. We can use this, then, to effectively iterate over the Cartesian product. Because the rtable is created with the storage=sparse option, the memory allocation is minimal. Here is a test procedure that implements this technique.

use_scanblock := proc(f)
local a,dims,n,L;
    L := args[2..-1]:
    n := nargs-1;
    dims := seq(1..nops(a), a in L);
    rtable_scanblock(rtable(dims, 'storage=sparse')
                     , [dims]
                     , proc(val,indx,data)
                       local i;
                           f([seq(L[i][indx[i]], i=1..n)]);
                       end proc
                    );
    NULL;
end proc:

use_scanblock(f,lsts): printf("\n"):
  [1, a]  [1, b]  [2, a]  [2, b]  [3, a]  [3, b]

Comparison

The following tabulates, for each procedure, the time and memory required to iterate through each element of a Cartesian product of lists. Four different inputs are tested. In each case, the rtable_scanblock technique is faster (about 2×) and uses less memory than combinat[cartprod].

Iterate through the Cartesian product of N lists of E elements,
generating E^N iterations.  At each step, call the (unassigned)
global function f with the generated list.

gcfreq = 10^7

               +-------- timing -------+--------- memory ---------+
256^2          | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     1.77  |     1.90  |    29.98   |     23.75   |    0
use_scanblock  |     0.93  |     0.95  |    19.23   |     13.87   |    0
---------------+-----------+-----------+------------+-------------+--------
10^5           |           |           |            |             |        
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     3.38  |     3.47  |    47.84   |     31.74   |    1
use_scanblock  |     1.84  |     1.93  |    28.49   |     22.68   |    0
---------------+-----------+-----------+------------+-------------+--------
3^10           |           |           |            |             |        
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     2.73  |     2.81  |    37.70   |     31.49   |    0
use_scanblock  |     1.40  |     1.45  |    21.51   |     16.06   |    0
---------------+-----------+-----------+------------+-------------+--------
2^16           |           |           |            |             |        
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     3.99  |     4.13  |    49.50   |     32.62   |    1
use_scanblock  |     2.06  |     2.08  |    25.63   |     21.00   |    0
---------------------------------------------------------------------------

Complete Instantiation

While iterators are useful for reducing the memory footprint, at times it is necessary to instantiate the Cartesian product as a single entity. Maple, alas, does not provide a function for doing so (this might make be a nice addition to the ListTools package). It is not difficult to create a procedure, the interesting part is making it fast. The following site has a few examples of implementations; I have used them below with minor modifications, along with some of my own contributions.

The first example uses combinat[cartprod], called Π_i n_i times, to generate the Cartesian product.

CartProdDeiterated := proc(L::seq(list))
local C,i,S;
    C := combinat['cartprod']([L]);
    [seq(C['nextvalue'](), i=1..mul(nops(S), S=[L]))]
end proc:

The next example uses a modification of the rtable_scanblock technique. Instead of calling rtable_scanblock, it generates an rtable with an initializer that maps each index to the corresponding entry in the Cartesian product. The elements of the rtable are then converted to a list by using seq.

CartRtable := proc(L::seq(list))
local i,a,n;
    n := nargs;
    [seq(a, a in rtable(seq(1..nops(a), a in args)
                        ,() -> [seq(L[i][args[i]],i=1..n)]
                        ,'order = C_order'
                       ))];
end proc:

Here is an offering from Carl DeVore (aka Carl Love). I've modified it slightly to use the seq parameter modifier previously mentioned. This routine produces the Cartesian product in column-major order; that is, the left-most element in the lists varies the fastest. All other routines produce the elements in row-major order.

CartProdAsList := proc(L::seq(list))
option `Copyright (C) 2002, Carl DeVore.  All rights reserved.`;
local q,i,N;
    N := map(nops, [L]);
    [seq([seq(L[i][1+irem(q,N[i],'q')], i= 1..nargs)], q=0..`*`(N[])-1)]
end proc:

Here is a suggestion from Helmut Kahovec, also slightly modified.

CartProdRecursive := proc(L::seq(list))
local i,j;
option `Copyright (C) 2001, Helmut Kahovec. All rights reserved.`;
    if nargs=2 then
        [ seq(seq([`if`(i::list,i[],i),
                   `if`(j::list,j[],j)
                  ] , j in args[2] )
              , i in L[1] ) ];

    elif nargs>2 then
        procname( procname(args[1..2]), args[3..-1] );
    else
        args;
    end if;
end proc:

The previous code has a few minor flaws and permits a few optimizations. Here is a non-recursive version that corrects the flaws and improves its performance.

CartProdNonrecursive := proc(L::seq(list))
local i,j,X,lst;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
    if nargs=0 then return [];
    elif nargs=1 then return map(`[]`,L);
    end if;
    X := [seq(seq([i,j], j in L[2] ), i in L[1] )] ;
    for lst in L[3..-1] do
        X := [seq(seq([i[],j], j in lst ), i in X )];
    end do;
    X;
end proc:

Here is an interesting version that I originally devised as a way to write an add-like procedure that summed over multiple dimensions, converted here to use seq to generate the Cartesian product. Because foldl is not a built-in, it is possible to improve it further, however, the resulting speed gains were trivial (less than five percent).

CartProdSeq := proc(L::seq(list))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
    eval([subs(Seq=seq, foldl(Seq
                              , [cat(i,1..nargs)]
                              , seq(cat(i,j)=L[j],j=nargs..1,-1)
                             ))]);
end proc:

To better understand what it does, run it without the eval statement, or do

test := subs(eval=''eval'',op(CartProdSeq)):
test([1,2,3],[a,b]); %;
            eval([seq(seq([i1, i2], i2 = [a, b]), i1 = [1, 2, 3])])
             [[1, a], [1, b], [2, a], [2, b], [3, a], [3, b]]

Comparison

The following tabulates, for each procedure, the time and memory required to generate the Cartesian product of a sequence of lists. Four different inputs are tested. For all cases the CartProdSeq is the fastest and uses the least memory.

Compute the cartesian product of N lists of E elements
to get a list of E^N sublists of N elements

gcfreq = 10^7

                     +-------- timing -------+--------- memory ---------+
256^2                | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     1.30  |     1.33  |    23.06   |     19.50   |    0
CartProdRtable       |     0.54  |     0.55  |    10.00   |      7.50   |    0
CartProdAsList       |     0.77  |     0.86  |    14.81   |     10.81   |    0
CartProdRecursive    |     0.25  |     0.26  |     7.79   |      5.75   |    0
CartProdNonrecursive |     0.13  |     0.13  |     5.77   |      3.75   |    0
CartProdSeq          |     0.12  |     0.13  |     5.76   |      3.75   |    0
---------------------+-----------+-----------+------------+-------------+--------
10^5                 |           |           |            |             |        
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     2.82  |     2.94  |    43.31   |     32.74   |    1
CartProdRtable       |     1.30  |     1.36  |    19.85   |     14.93   |    0
CartProdAsList       |     1.94  |     1.98  |    27.67   |     21.93   |    0
CartProdRecursive    |     0.74  |     0.74  |    18.41   |     13.56   |    0
CartProdNonrecursive |     0.47  |     0.49  |    14.78   |      9.94   |    0
CartProdSeq          |     0.32  |     0.35  |    13.53   |      9.37   |    0
---------------------+-----------+-----------+------------+-------------+--------
3^10                 |           |           |            |             |        
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     2.28  |     2.34  |    31.04   |     27.99   |    0
CartProdRtable       |     1.08  |     1.10  |    12.54   |     10.56   |    0
CartProdAsList       |     1.85  |     1.87  |    24.33   |     20.68   |    0
CartProdRecursive    |     0.66  |     0.73  |    18.49   |     13.37   |    0
CartProdNonrecursive |     0.43  |     0.46  |    14.54   |     10.19   |    0
CartProdSeq          |     0.24  |     0.27  |    10.56   |      8.25   |    0
---------------------+-----------+-----------+------------+-------------+--------
2^16                 |           |           |            |             |        
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     3.67  |     3.80  |    44.03   |     35.24   |    1
CartProdRtable       |     1.64  |     1.67  |    17.02   |     14.56   |    0
CartProdAsList       |     3.06  |     3.18  |    35.65   |     30.74   |    0
CartProdRecursive    |     1.15  |     1.31  |    28.41   |     23.31   |    0
CartProdNonrecursive |     0.78  |     0.83  |    23.88   |     19.25   |    0
CartProdSeq          |     0.44  |     0.48  |    18.21   |     14.00   |    0
---------------------------------------------------------------------------------

Frontend

The frontend command is a bit tricky. The basic idea is that it replaces expressions in the arguments to a procedure with names, evaluates the procedure with the replaced arguments, then back-substitutes the expressions for the the names in the final result. The first argument passed to frontend is the procedure that is to be evaluated. The second argument is a list containing the arguments to be passed to the procedure. Remaining arguments are optional, we'll get to those in a second. Suppose that we are given the expression

y := cos(x) + 3*a*x^2:

and we want the coefficient of the Maple Equation term . We might try using coeff:

coeff(y, x^2);

Error, unable to compute coeff

That fails because y is not a polynomial in x. The frontend command can be used here:

frontend(coeff, [y,x^2]);

To get a better handle on how this works, I'll create a modified version of frontend, called Frontend:

Frontend := proc(f)
frontend( proc()
printf("%a\n",'f(args)');
f(args) end proc, _rest)
end proc:

This version prints what the actual call looks like, before evaluating it:

Frontend(coeff, [y,x^2]);

coeff(O+3*a*x^2,x^2)

Here we see that the actual call to coeff is coeff(O+3*a*x^2,x^2). This is clearly a polynomial in x, so coeff is able to return the desired result. You are probably wondering what the 'O' is doing in the expression. That is a frozen version of the expression cos(x). Let's try a slightly different example:

Frontend(coeff, [cos(x)+sin(x)+5*x^2,x^2]);

coeff(O+O+5*x^2,x^2)

This time the first argument passed to coeff is O+O+5*x^2. The two O's are frozen expressions for cos(x) and sin(x). While they look the same, they are actually different names to Maple.

Now let's examine the optional third argument to frontend. It must be a list of exactly two sets. By default, this argument has the value [`+`,`*`,], which means that sums and products are not frozen. The first set in the list is a set of types. Any subexpressions in the list of arguments that match these types are not frozen. Let's see this in action:

Frontend(coeff, [a*cos(x)+b*sin(x)+c*x^2, x^2], [`+`,]);

coeff(O+O+O,x^2)

Unlike the default, subexpressions of type `*` are frozen. That leads to an error, since the c*x^2 term was frozen. One way to avoid that is to include the term x^2 in the second set of the third argument. Expressions in the second set are not frozen (the first set consists of types not to freeze, the second set consists of expressions not to freeze.

Frontend(coeff, [a*cos(x)+b*sin(x)+c*x^2, x^2], [`+`,x^2]);

coeff(O+O+c*x^2,x^2)

Now try calling frontend with the coefficient of x^2 set to 1 and without including x^2 in the set of expression not to freeze.

Frontend(coeff, [a*cos(x)+b*sin(x)+x^2, x^2], [`+`,]);

coeff(O+O+x^2,x^2)

Notice that the x^2 term was not frozen. frontend never freezes, as a single expression, an exponent of an integer. That includes an expression of the form 1/expr, which is represented internally in Maple as expr^(-1).

With that in mind, we now return to the original problem, which was computing the eigenvectors of a Matrix of functional expressions. First, let's pass use the default third argument:

Frontend(LinearAlgebra:-Eigenvectors,[<<f(a)|0>,<0|f(b)>>]);

LinearAlgebra:-Eigenvectors(O)

Error, (in LinearAlgebra:-Eigenvectors) invalid input: LinearAlgebra:-Eigenvectors expects its 1st argument, A, to be of type Matrix(square) but received O

That didn't work. We need to prevent frontend from freezing the complete Matrix. To do that, include the type Matrix in the set of types that are not to be frozen.

Frontend(LinearAlgebra:-Eigenvectors,[<<f(a)|0>,<0|f(b)>>],[Matrix,]);

LinearAlgebra:-Eigenvectors(Matrix(2, 2, [[O,0],[0,O]]))

That works. Consider a slightly more complicated case:

Frontend(LinearAlgebra:-Eigenvalues,[<<f(a)+f(b)|f(b)>,<f(a)|f(b)>>],[Matrix,]);

LinearAlgebra:-Eigenvalues(Matrix(2, 2, [[O,O],[O,O]]))

That is correct, but the result is more complicated than necessary because frontend froze the polynomials inside the Matrix. To prevent that we can include the sum and product types in the set of types not to freeze, along with Matrix:

Frontend(LinearAlgebra:-Eigenvalues,[<<f(a)+f(b)|f(b)>,<f(a)|f(b)>>],[Matrix,`+`,`*`,]);

LinearAlgebra:-Eigenvalues(Matrix(2, 2, [[O+O,O],[O,O]]))

Maple Equation

This post was generated using the MaplePrimes File Manager

View 1_Frontend.mw on MapleNet or Download 1_Frontend.mw
View file details

Functional Object Programming with Maple

As you can see, this technique implements quite full and simple functional objects for Maple.

It is very simple and efficient. All properties are private, you need to implement public get and set methods to make them public.

Public methods and stored in the last table structure of the procedure.

I'm deeply interested by comments...

retart;

Maple Equation

Functional object programming with Maple

Simplest functional object

simplest_fo:=proc(n) local this, set_;
this := [n];
set_:=proc( p); this[1]:=p; end proc;
table(['get'= (x->this[1]), 'set'= set_]);
end proc;

Maple Equation

foo:=simplest_fo( 2);

Maple Equation

foo['get']();

Maple Equation

foo['set'](5);

Maple Equation

foo['get']();

Maple Equation

Application to probability distributions

G:=(m,s)->(t->(exp(-(t-m)^2/(2*s^2))/(sqrt(2*Pi)*s)));

Maple Equation

structured_distribution:=proc( f, A, B)
table(['density'= (x->f(x)), 'repartition'=(x->Int(f(t),t=A..x)), 'plot'=(t->plot(f(x),x=A..B))]);
end proc:

SD:=structured_distribution( G(0,1), -infinity, infinity):

SD['density'](1);

Maple Equation

SD['repartition'](1);

Maple Equation

SD['plot']();

Maple Plot

This post was generated using the MaplePrimes File Manager

View 744_indexed-output.mw on MapleNet or Download 744_indexed-output.mw
View file details

Maple 401

So you have used Maple as a glorified calculator (Maple 101), then wrote a few 1 liners (Maple 201), and even a few larger procedures (Maple 301), where you were both amazed and horrified by 'op'. But when you get serious about programming in Maple, even for not-so-large procedures, what are the fundamental parts of the system that you should know?

Other pages in this book talk about particular features. This one is instead a simple list of those Maple commands and concepts you need to know to be able to call yourself a Master Maple Programmer.

Concepts you need to understand:

structured types (see ?type,structured)
DAGs versus trees
hash tables
expression sequences (and why they are both useful and dangerous)
rtables and indexing functions
evaluation rules (one level, last-name-eval, special rules)
higher-order functions
lexical scoping
environment variables (known as fluid variables in LISP)
argument processing
attributes
why you should never use 'whattype'
why you should be wary of 'with'
why 'has' and 'hastype' can kill efficiency
options Cache, remember, system, builtin, inlne, etc
what all the quotes mean
the difference between evaluation and automatic simplification

Commands you need to know:

indets
map, zip and seq
member, select, remove, selectremove
`if`
unapply
the 'in' operator
2-argument eval
procname, args, nargs, and their friends
frontend, freeze, thaw
maptype
depends
use, uses and bind
kernelopts and interface
typematch
subsop
unames, anames
LibraryTools, march
$define, $include, $ifdef and other pre-processor macros
system, ssystem, process
evalindets and subsindets
ToInert and FromInert
call_external
dismantle, disassemble, pointto, addressof

Measuring Maple with Macros

This book entry describes a method, using text files and preprocessor macros, to measure the timing and memory usage of Maple procedures. Attached are the source code and compiled library for a small Maple package that does the actual measurements.

Introduction

One of the goals, possibly the main goal, of this web book, High Performance Maple Programming Techniques, is to demonstrate methods for increasing the speed of Maple computations. To determine whether a particular algorithm is faster than another, the execution time of the two must be accurately measured. Speed, however, is not the sole criterion for judging performance. Memory usage is also a significant metric, particularly if the algorithm must work with large data sets. This chapter describes a method for measuring the execution times and memory usage of competing Maple procedures.

The primary difficulty in accurately measuring the performance of several procedures is ensuring that the Maple enviroment remains constant for each test run. The best way to achieve that is to start a new Maple session for each measurement; that, however, makes writing the test code and comparing the results a nuisance.

The next best approach is to issue a restart command between each test run. That causes Maple to clear its internal memory and return to the operating system most of the memory that it has allocated. This can be done in the Maple worksheet environment, however, because the code required to initialize the data must be duplicated for each procedure that is tested, this method can be tedious and unwieldy. For example, changing the initial data requires editing the same code in several places in the worksheet. Because each restart command clears internal memory, collecting the results into a single table so that they can be quickly compared requires writing the results of each test to a file and then reading them back into the worksheet. Finally, while the worksheet enviroment is excellent for interacting with Maple and for presenting Maple results, a programmer's text editor (emacs, vi, jedit, etc) is better for writing Maple code.

The following sections describe a method that uses tty maple (command-line maple) to generate a table comparing the execution time and memory usage of alternative procedures. All code being tested may be stored in a single text file. The initialization code that is executed need only be entered in one place, so it can be readily modified. To accomplish this, Maple preprocessor macros are used to recreate the initialization data between restart commands. The use of Maple preprocessor macros precludes using the Maple worksheet interface for the tests.

Example

Test File

The code below shows part of the test file used to measure the performance of procedures for sorting list of names. The complete file is called sort-names.mpl and is attached to this web page. Explanations of each section follow.

      #-------- Define Preprocesser Macros ----------
      $define LOOPS 200 # number of time procedure is executed
      $define NUM 50    # number of words in list
      $define LENGTH 5  # number of characters in each word
      $define ARGS [seq](StringTools:-Random(LENGTH, 'alpha'), cnt=1..NUM)

      $define MEAS(i) \
        kernelopts('gcfreq'=10^7): \
        StringTools:-Randomize(1): \
        Usage:-TimeMemoryPrint(SortNames[i], argseq = ARGS, loops = LOOPS ): \
        restart;

      #-------- Print Notes and Table Header -----------
      printf("Testfile: %s\n\n", TESTFILE):
      printf("Compare performance of sorting names.\n"):
      printf("Each test sorts the same list of %d names of %d characters %d times.\n\n", NUM, LENGTH, LOOPS):
      Usage:-TimeMemoryHeader();

      restart;

      #------- Assign and Measure each Procedure -------

      SortNames[2] := proc(names :: list)
          sort( names
                , (a,b) -> sprintf("%a",a) <= sprintf("%a",b)
              );
      end proc:

      MEAS(2);

      SortNames[4] := proc(names :: list)
      local n;
          map(attributes, sort([seq](setattribute([sprintf]("%a",n),n)
                                     , n = names)
                               , 'lexorder'[1]));
      end proc:

      MEAS(4);

      done

When executed with tty maple this produces the following output:

      Testfile: sort-names.mpl

      Compare performance of sorting names.
      Each test sorts the same list of 50 names of 5 characters 200 times.

      time (s)  used (MB)  alloc (MB)  gctimes  proc
      --------  ---------  ----------  -------  ----
        1.43     18.89       14.37         1    SortNames[2]
        0.17      2.51        1.50         1    SortNames[4]

Macro Definitions

Maple preprocessor macros are used to control the execution of the tests. Macros are used rather than Maple variables because these must survive a Maple restart command.

The LOOPS, NUM, and LENGTH macros correspond, respectively, to the number of repetitions the test procedure is executed, the number of names in each list, and the length of each name.

The ARGS macro generates the list of names, using the defined macros for parameters. The resulting list is regenerated and passed to each test procedure during measurements.

The MEAS macro has one parameter, the index of the test. The macro expands to a series of Maple statements that are executed for each measurement. The following describes the effect and purpose of each statement.

The call to kernelopts('gcfreq'=10^7) reduces the occurrence of garbage collection.
StringTools:-Randomize(1) ensures that the same list of names is generated for each test.
Usage:-TimeMemoryPrint(...) passes the expression sequence contained in the list [ARGS] to SortNames[i] and executes it LOOPS times. It computes the time and memory usage and prints a single row in the output table.
The restart command restarts Maple. This ensures that each test starts with the same environment.

Notes and Table Header

The first printf statement makes use of the preprocessor macro TESTFILE, which is assigned by tmaple (see below) to print the name of the file being tested The next two printf statements print the common conditions for the test. The call to Usage:-TimeMemoryHeader() prints the table header and the row of underscores beneath it. A restart command before the next part ensures that each test seems almost the same starting condtions.

Procedure Assignments and Measurements

Each procedure to be measured (Sort[2] and Sort[4]) is assigned and then measured using the previously defined macro MEAS.

Usage Package

The Usage package is a small, custom package that exports commands for measuring and tabulating the time and memory usage of procedures being tested. The source code, Usage.mpl, and compiled Maple archive, Usage.mla, are attached at the bottom of this web page. The package exports three procedures:

TimeMemoryGC(ex :: uneval, loops :: posint := 1)
Evaluates the expression ex in a loop loops times and returns a sequence of four expressions:
1. the duration (in seconds) of the loops to evaluate
2. the additional memory used (in bytes)
3. the additional memory allocated (in bytes)
4. the number of garbage collections.
TimeMemoryPrint(p, {argseq := NULL, loops :: posint := 1} )
Calls TimeMemoryGC( p(argseq), loops ) and prints the result as a row of a table. The arguments to the procedure and the number of loops are keyword parameters (a new feature as of Maple 10); their defaults are NULL and 1, respectively. It returns NULL.
TimeMemoryHeader
Prints a header for the table. It returns NULL.

Maple Script

After installing the Usage.mla archive in a directory that is in your libname path, you can measure the performance of, for example, the test file sort-names.mpl by issuing the command

      $ maple -q sort-names.mpl

Windows users should use the cmaple command, see the Maple help page ?maple for details.

While that works, a possible disadvantage to calling tty maple directly is that it normally loads the user's Maple initialization file, which may assign additional, unneeded procedures. A convenient alternative is to execute a shell script that calls tty maple using specific options to avoid reading the initialization file. The attached file tmaple does just that. Note that this is a Linux shell script; it should be usable in Windows via cygwin.

To use the shell script, install it in your path and do

      $ tmaple sort-names.mpl

Calling it with the -h option prints a short list of available options.

Attachments

tmaple.txt This file should be renamed tmaple, the txt extension was necessary to upload it to the MaplePrimes site.
sort-names.mpl Sample test file.
Usage.mpl Source code for Usage.mla.
Usage.mla Maple archive.

Sorting Lists of Names

This book page describes and compares several methods for sorting lists of names. It revisits using attributes for rapid sorting, with attention to their global effect on strings. It demonstrates the use of the lexorder[n] option to sort and it illustrates the use of convert/local.

Introduction

Occasionally it is useful to sort a list of names—usually variables extracted from an expression—into a repeatable order. The simple and obvious way to do this is with the sort command.

      sort([x, y, Z, a1, a1]);
                                     [Z, a1, a1, x, y]

That approach does not work if there are indexed names in the list.

      sort([x, y, Z, a[1], a[2]]);
                                    [x, y, a[1], a[2], Z]

The elements are sorted, but by machine address rather than alphabetically. This can be shown using the addressof command:

      map(addressof,%);
                [135350720, 141574716, 143418428, 143418444, 143451544]

The explanation can be found from a careful reading of the sort help page. A list of symbols are, by default, sorted lexically. However, an indexed name is not a symbol. Consequently, the list is sorted by machine address, which is session dependent and hence not repeatable.

Sorting Techniques

The following sections look at a few ways to use the Maple sort command to sort a list of names into a repeatable order. Ideally the same final permutation is returned for any given permutation of the elements in the list.

Sorting as strings

The sort command takes an optional second argument which can be an ordering function that takes two arguments—elements of the list—and returns true if the first argument is sorted with respect to the second, false otherwise. To sort names we can compare them as strings.

      SortNames[1] := proc(names :: list)
          sort( names
                , (a,b) -> convert(a,string) <= convert(b,string)
              );
      end proc:

      SortNames[1]( [x, y, Z, a[1], a[2]] );
                                             [ Z, a[1], a[2], x, y ]

Note the use of the inequality `<=' in the comparison rather than strict inequality, `<'. That prevents an unnecessary swap if the strings are equal, and achieves a stable sorting, that is, one in which equal elements maintain their original order.

This procedure has two problems. First, it is rather slow, at least in comparison to later techniques. Second, because they map to the same string, it cannot distinguish the indexed variable x[1] and the symbol `x[1]`. We can correct this latter problem by using sprintf with the %a format string to convert names to strings:

      SortNames[2] := proc(names :: list)
          sort( names
                , (a,b) -> sprintf("%a",a) <= sprintf("%a",b)
              );
      end proc:

      SortNames[2]( [x, y, Z, `a[1]`, a[1], a[2], `b[1]` ] );
                                    [Z, `a[1]`, `b[1]`, a[1], a[2], x, y]

Sorting with attributes

In a previous blog entry

I explained how attributes can be used to sort a list of floating-point numbers with respect to a relatively expensive function. That speeded up the operation by (1) reducing the number of calls to the function to the number of elements in the list, and (2) removing the need to call a user-defined order function, which is invariably slower than the built-in Maple order functions. The same technique can be employed here. However, because we are sorting strings rather than floats, there are a few new subtleties to consider.

Here is the basic routine

      SortNames[3] := proc(names :: list)
      local n;
          map(attributes, sort([seq](setattribute(sprintf("%a",n),n)
                                     , n = names)
                               , 'string'));
      end proc:

If you have not read the previously mentioned article this may need a bit of explanation. The [seq](...) expression creates a list of strings, using sprintf as previously described. As each string is created it is passed to setattribute which assigns the original name as the attribute to the string. The resulting list of strings is then sorted lexically. Finally the attributes of the sorted strings are extracted and returned. Because the attributes are the original names of the strings, the returned list is the sorted list of names.

This procedure has two flaws: one minor, one major.

The minor flaw is that it has global side-effects due to the way that Maple handles strings. Each unique string in Maple exists in memory in one location. Adding an attribute to a string does not change this location. Calling setattribute replaces any preexisting attributes, so if some other procedure was using attributes with strings, and one of the strings matched a name passed to this procedure, its attribute would be changed.

The major flaw is subtle but potentially fatal. Suppose two distinct names mapped to the same string. When the attributes are extracted from the sorted strings, both strings are mapped to the same name, the second one in the original list. The returned list of names, then, does not contain one of the original names. Using sprintf rather than convert/string eliminates some of the collisions, however, it is still possible for two names to map to the same string. A simple example is escaped local variables, which can be readily generated with convert/local:

      lst := map(convert, [a,a], `local`):
      convert(lst, set);
                                                 {a,a}

To demonstrate that this really occurs convert the sorted list to a set.

      convert(Sort[3](lst), set);
                                                  {a}

Avoiding Collisions

Both flaws in SortNames[3] can be corrected by attaching the attributes to a list containing a string rather than the string itself. Unlike attributes attached to a string, attributes attached to a list (and most attributable objects) only affect that instance of the object; two lists that are otherwise identical but have separate attributes are different elements in Maple. This can be verified with

      lst := [1,2]:
      lst1 := setattribute(lst,1):
      lst2 := setattribute(lst,2):
      lst3 := setattribute(lst,2):
      map(evalb, [lst1=lst2, lst2=lst3]);
                                          [ false, true ]

Normally, enclosing elements of a list in lists requires providing a user-defined order function to sort to unpack and compare them, which reduces the advantage of this technique. However, sort takes the optional symbol lexorder[n] which tells it to lexically compare the n'th element of each sublist. The following procedure does just that, using n=1.

      SortNames[4] := proc(names :: list)
      local n;
         map(attributes, sort([seq](setattribute([sprintf]("%a",n),n)
                                    , n = names)
                              , 'lexorder'[1]));
      end proc:

Demonstrate that this fixes the collision problem.

      convert(SortNames[4](lst), set);
                                       {a,a}

Sorting as Symbols

The procedures discussed convert names to strings and then sort the strings. Another option is to convert names to symbols and sort the symbols. The Maple procedure convert/symbol may be used to convert names (and other expressions) to symbols. However, the distinct names a[1] and `a[1]` get mapped to the same symbol. To avoid that, we can use convert/local, which maps each name to a distinct symbol.

      SortNames[5] := proc(names :: list)
      local n;
          map(attributes, sort([seq](setattribute(convert(n,'`local`'),n)
                                     , n = names)));
      end proc:

      convert(SortNames[5](lst), set);
                                       {a,a}

Remaining Issue

The goal was to create a procedure that rapidly sorts any permutation of a list of names into a session-independent, repeatable, sorted order. The results, SortNames[4] and SortNames[5], do not quite meet that goal. If the list of names contains names that lexically match but are different (escaped local variables) then the resulting order of those names is identical to their relative input order. Otherwise those procedures work as required. This is a relatively benign failure; furthermore it is not clear it can be be avoided.

Performance Comparison

As shown in the following table, using attributes (cases 3-5) is substantially faster and uses less memory than providing a user-defined order function to sort (cases 1-2). SortNames[4], which uses lists to avoid string collisions, is slighter slower, but more robust, than SortNames[3], but faster than SortNames[5], which uses symbols. SortNames[4] is the best overall choice.

      $ tmaple sort-names.mpl

	Compare performance of sorting names.
	Each test sorts the same list of 50 names of 5 characters 500 times.

	time (s)  used (MB)  alloc (MB)  gctimes  proc
	--------  ---------  ----------  -------  ----
	  2.28     54.70       36.12         2    SortNames[1]
	  3.73     47.19       35.93         2    SortNames[2]
	  0.40      4.90        3.44         1    SortNames[3]
	  0.46      6.25        3.50         1    SortNames[4]
	  0.60     10.63        6.56         1    SortNames[5]

The table was generated using custom tools for accurately and repeatably measuring the time and memory usage of Maple procedures. They shall be the subject of another article.

Sparse Matrix-Vector product using evaluation

This tip comes care of Dr. Michael Monagan at Simon Fraser University. Represent your sparse matrix as a list of rows, and represent each row as a linear equation in an indexed name. For example:

A := [[1,0,3],[2,0,0],[0,4,5]];
S := [ 1*x[1] + 3*x[3], 2*x[1], 4*x[2]+5*x[3] ];

To compute the product of the matrix A with a Vector X, assign x[i] := V[i] and evaluate. This can be done inside of a procedure because x is a table.

V := [7,8,9]:
for i to 3 do x[i] := V[i] end do:
eval(S); # result
for i to 3 do x[i] := evaln(x[i]) end do: # unassign

Don't forget to unassign. Note that this method of assignment is slow when V is a big list. In short, we need to use "for i in V", which does one nice linear pass through V. This is what I do:

n := 0;
count := proc() global n; n := n+1 end proc:
for i in V do x[count()] := i end do:
eval(S);
for i to 3 do x[i] := evaln(x[i]) od:

In a real procedure, n would be lexically scoped and not a global variable. Instead of an explicit loop you can also use seq and assign.

n := 0;
count := proc() global n; n := n+1 end proc:
seq(assign(x[count()],i), i=V);
eval(S);
seq(unassign('x[i]'), i=1..3);

Of course, where this trick really shines is in doing back substitution. You can easily substitute tens of thousands of variables into any number of equations, all in linear time.

The Top Ten Maple Errors

I've made up a worksheet of the Top Ten Maple Errors, containing some of the common mistakes I often see newcomers to Maple commit (especially in the setting of my Introduction to Mathematical Computing class). I hope you will find it useful in trying to avoid those mistakes. Of course this is only a personal list, and not exhaustive. Please feel free to argue the merits of other items that should be included in the list.

Here is the link:

View 4541_topten.mw on MapleNet or Download 4541_topten.mw
View file details

Tine: A Module for Measuring Maple Performance

Tine: A Module for Measuring Maple Performance A previous entry in this book, Tools of the Maple Masters, described an idiosyncratic technique for measuring the time and space (memory) usage incurred while evaluating expressions in Maple. It used Maple preprocessor macros so that the restart command could be executed between measurements. As explained there, restarts are used to ensure that each evaluation starts with a minimal usage of memory. While that method was suitable for measuring a procedure at a few data points, the lack of looping with the available macros made it impractical for more intensive analyses. Furthermore, macros are only usable with tty Maple (cmaple), are painful to debug, and usually involve creating two files: one that assigns the procedures to be tested, the other that runs the tests.

I wanted something better. In particular, I wanted to be able to easily plot the performance of a procedure, or several procedures, as the size of the input increased, with just a few lines of Maple code. The key to doing this is controlling, and possibly restarting, an independent Maple process. Preliminary attempts, using the Sockets package and OpenMaple, to communicate with a Maple process, were not entirely successful. Each has limitations.

Communicating with an external Maple process is not difficult, the hard part is restarting it. The Sockets package, which uses network sockets

to communicate, does not provide a way to programmatically restart the kernel. Furthermore, once one end of a socket disconnects (for example, after killing the process) it is not possible, apparently, to reconnect. OpenMaple, which permits accessing Maple through compiled code, has a RestartMaple command. Alas, OpenMaple cannot be directly used by another Maple kernel, at least not via the ExternalCalling procedures. From the help page:

These functions [StartMaple, StopMaple, RestartMaple] can be used in external code with OpenMaple. They do not work as part of define_external.

During an evening walk, while considering methods for bypassing these obstacles, I realized a simple technique for accomplishing the task. I had taken the left fork of a path that led up a hill back to my home, noticed that it was a bit steeper than I cared to climb, so chose an alternate path that returned to the main fork. The technique I wanted was analogous to my path, it uses a process fork.

Caveats

Before describing the method and implementation, there are a few caveats to be considered. First, forks are available only on some operating systems. They do not work on Windows systems. Second, as of Maple 11, the process package is deprecated. From the process[fork] help page:

The process package has been deprecated for Maple 11. The fork command may be removed from future versions of Maple. For similar functionality see the Threads package.

Alas, I do not believe that the Threads package provides equivalent functionality. A third caveat is that, to increase the utility, you will need to compile a small external C program. A makefile is provided, but I have only used it on my Linux system. The package works without the external program, but having it increases the precision of the timing measurements.

Description

A Unix process fork causes the current process to fork into two branches: the parent branch, which we identify with the original process, and the child branch, which we consider a new process. At the time of the fork, the two processes are identical, but independent. Well, not completely identical. A variable is set to one value (0) in the child process and to another value (the process ID, or pid, of the child) in the parent process. The code in each process inspects this variable and branches accordingly.

To use a fork to measure the performance of a Maple procedure, we first assign the procedure and variables in the main branch, open a pipe, then fork the process. The child process records the time and memory usage, applies the test procedure to the arguments, measures the time and memory usage, and computes the differences. It converts the measured data to a string, transmits it on the pipe, then kills itself.

Meanwhile, the parent process reads the other end of the pipe. When the transmitted data arrives, it closes the pipe, waits for the child process to terminate, parses the string, and returns the measured data as the output. Because any memory allocated by Maple during the evaluation of the test procedure only affects the child process, the memory usage in the parent process is minimal—this ensures that the child process started by the next measurement has essentially the same start conditions. Killing the child process also permits the OS to reclaim its memory, something that is not done when Maple is merely restarted.

Implementation

The Maple module Tine exports several procedures:

ApplyFork: Apply a procedure to arguments in a child process and return the result.
ApplyMeas: Measure the time, memory usage, and number of garbage collection cycles required to apply a procedure to arguments.
EvalFork: Evaluate an expression in a child process and return the result.
EvalMeas: Measure the time, memory usage, and number of garbage collection cycles required to evaluate an expression.
PlotApplyMeas: Plot the measurements of a procedure, or procedures, as a parameter varies.
SeqLog: Return a logarithmic sequence. Useful for generating domain points that are evenly spaced on a logarithmic scale.
Time: Return the time, with microsecond resolution, since the epoch. This is available only if the accompanying C code has been properly compiled.

Example

Compare catenating strings with the cat function, the || operator, and StringTools[Join]. The use_catbin function uses cat as a purely binary procedures, to give a useful comparison with ||.

# Assign procedures to measure.  Each procedure catenates all its arguments.

use_bar    := proc() foldl(`||`,_passed) end proc:
use_catbin := proc() foldl(cat,_passed) end proc:
use_join   := proc() StringTools:-Join([_passed],"") end proc:
use_cat    := proc() cat(_passed) end proc:

# Load the Tine package
with(Tine):

# Assign N a list of integers, [n_i], logarithmically spaced.
# These are used to generates sequences of that many strings.
# The Loops variable is assigned a list of integers approximately
# inversely proportional to the corresponding integers in N; they are
# used to set the number of times the test procedure is evaluated.
# That permits reasonably accurate measurements for small n_i.

(N,Loops) := SeqLog(10..10^3,30,andinverse):

# Generate a list of lists of strings to catenate.

Pargs := [seq([cat("a",1..n)], n=N)]:

# Measure and plot the bytes-used and total-time required by each of the
# procedures.  The 'seed = true' argument causes an integer to be appended
# to each sequence of arguments passed to the test procedures so that the
# results from consecutive calls, when multiply evaluated, differ.
 
PlotApplyMeas([use_bar, use_catbin, use_join, use_cat]
              , N, Loops, Pargs, 'seed=true'
              , field = [bytesused, totaltime]
             );

Installation

The package is distributed as a zip file. It contains a README file that briefly explains how to install it.

Turn loops into seqs with assign

Consider the following three problems:
1) given a list [a,b,c,d,e], return the list [a=1,b=2,c=3,d=4,e=5]
2) given a nested list [a, [b, [c, [d, [e]]]]] return the list [a,b,c,d,e]
3) given an integer in base 10, compute it's base b representation

First I will show you what not to do:

L := [a,b,c,d,e];
M := [];
for i from 1 to nops(L) do
  M := [op(M), L[i]=i];
end do;

Building up lists (and sets) incrementally is quadratic time, because each iteration of the loop allocates linear storage to hold the new list. The standard solution is a loop with a temporary variable, assigning to a table:

M := table():
for i to nops(L) do
  M[i] := L[i]=i;
end do:
[seq(M[i], i=1..nops(L))];

However for large lists this is slow, because you are accessing the elements of L incrementally.
If you time it, you will see that seq(L[i], i=1..nops(L)) is much slower than seq(i, i=L) when L is a large list. Thus we need to do one linear pass through the list L. We can modify the loop as follows to acheive the desired result:

M := table():
j := 1:
for i in L do 
  M[j] := i=j;
  j := j+1;
end do:
[seq(M[i], i=1..nops(L))];

The code above is asymptotically faster than all previous versions, however it still requires a table and random accesses to construct the result. Can't we just somehow use seq ?

The answer is "yes", and the trick is to use the assign command: assign('x', v) assigns v to x and returns NULL. Here is the loop above, using one seq and two local variables i and j:

L := [a,b,c,d,e];
j := 1;
[seq(i=`+`(j, assign('j', j+1)), i=L)];

Notice that I had to add the NULL to j in order to make it all fit in a seq. If this fails, for example if j is a modp1 polynomial, then you can do op([j, assign('j', j+1)]) instead which is slower but more general.

Now lets look at the second problem. Nested lists are interesting because they are one way to implement stacks in Maple. The underlying lists [a, [...]] store a pointer to an object and a pointer to a substack. I don't know whether you actually want to use nested lists for stacks (you certainly don't want to display them), but the problems we face here pop up in other settings, for example you may need to extract linked list data from external code. Again the standard solution is a loop and table:

L := [a, [b, [c, [d, [e]]]]];
M := table():
i := 1:
for i while nops(L) > 1 do
  M[i] := L[1];
  L := L[2];
end do:
M[i] := L[1]:
[seq(M[j], j=1..i)];

The point is that after we extract the leading element we assign the sublist to L. Here it is with a seq:

L := [a, [b, [c, [d, [e]]]]];
L := [seq(`+`(L[1], assign('L', L[2])), i=1..4), L[1]];

Of course we cheated - how does one know ahead of time how many elements are in L ? The range for seq must be pre-specified and can not be changed, so the only solution for now is to guess and overshoot. Perhaps the Maple kernel developers will consider adding a seqwhile command in the future :) We will use an `if` statement to detect the end of the stack and return NULL.

L := [a, [b, [c, [d, [e]]]]];
L := [seq(`if`(L=[], NULL, `+`(L[1], `if`(nops(L)=1, assign('L',[]), assign('L', L[2])))), i=1..10)];

This does too many comparisons and is slow. You can speed it up by using an explicit delimiter, where the last list element is [e, []]:

L := [a, [b, [c, [d, [e, []]]]]];
L := [seq(`if`(L=[], NULL, `+`(L[1], assign('L', L[2]))), i=1..10)];

While we're on the topic, here's how to build that list with a seq. Note that the order must be reversed because we are building a stack.

L := [e,d,c,b,a];  # input reversed
j := []:
seq(assign('j', [i, j]), i=L):  # outputs nothing
j;  # result

Now for the final problem, which was the original motivation. I want to compute the base b representation of an integer with the sign as the first element. Normally one would use `convert/base` but in my case it was too slow. I had thousands of integers to convert. My solution was based on the following:

b := 11:
L := [4387324, -9879862, 4237654, -123489321];
# bound the number of base b digits required
k := trunc(evalf(ln(max(seq(abs(i), i=L)))/ln(b)))+1;
[seq([sign(i), assign('t',abs(i)), seq(irem(t, b, 't'), j=1..k)], i=L)];

# here is the code I replaced:
[seq([sign(i), op(convert(abs(i), base, b))], i=L)];

Note that the ceil function is very slow and should be avoided, along with floor, frac, and round. The irem trick is used by `convert/base` as well, however if you time the two methods on large lists you should find that my method is faster because there is much less overhead. Also, if your ultimate goal is to get a matrix of integers, you can use the (faster) rtable constructor instead of the Matrix command because the sublists are all the same size. In my case I was even able to compute the bound k must faster (using ilog2 and iquo) because my base was 2^32.

There are many other applications of this "assign in a sequence" technique, and I would be interested to hear what people might have tried.

Using try/finally to temporarily modify a global flag

Here is the article. I wrote this as a blog entry, and felt it was better to leave it one place rather than duplicating it. Besides, I don't know how to delete a blog entry. A cheap way to earn five points; however, my next page addition—in the works—to this book will be recompense.