Joe Riel

9660 Reputation

23 Badges

20 years, 13 days

MaplePrimes Activity


These are replies submitted by Joe Riel

Interesting:

    |\^/|     Maple 13 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2009
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> proc(x); local a; end proc;
                                                 proc(x) local a;  end proc

Your penultimate example is surprising; it indicates that the Standard 1D parser differs from the command-line parser.  That is

(**) proc(x); local a; end proc;

is syntactically acceptable in command-line Maple.

That is a good question. It works with cos because cos is a Maple procedure.  So you can just pass the name of the procedure to bisection. Similarly, you could have done

f := x -> x^2 - 2: # create a procedure
bisection2(f, 1, 3, 1e-7);

 

That is a good question. It works with cos because cos is a Maple procedure.  So you can just pass the name of the procedure to bisection. Similarly, you could have done

f := x -> x^2 - 2: # create a procedure
bisection2(f, 1, 3, 1e-7);

 

Did you try evaluating f(c) inside the debugger:

DBG> f(c);
x(2.000000000)^2-2

That should tell you something. What is happening is that Maple is applying an argument, 2.0, to the expression x^2-2 and computing (correctly, if unexpectedly) x(2.0)^2-2.  The problem is that you passed an expression rather than a procedure as the first argument to bisection2. To catch that error early on, one can declare f to be a procedure:

bisection2 := proc(f::procedure, A::realcons, B::realcons, delta::positive)
local a,b,c;
    a := evalf(A);
    b := evalf(B);
    while abs(b-a)>delta do
        c := (a+b)/2:
        if sign(f(c))=sign(f(a)) then
            a := c;
        else
            b := c;
        end if;
    end do;
    (a+b)/2;
end proc:

You should have done

bisection2(x -> x^2-2, 1, 3, 1e-7);

 

Did you try evaluating f(c) inside the debugger:

DBG> f(c);
x(2.000000000)^2-2

That should tell you something. What is happening is that Maple is applying an argument, 2.0, to the expression x^2-2 and computing (correctly, if unexpectedly) x(2.0)^2-2.  The problem is that you passed an expression rather than a procedure as the first argument to bisection2. To catch that error early on, one can declare f to be a procedure:

bisection2 := proc(f::procedure, A::realcons, B::realcons, delta::positive)
local a,b,c;
    a := evalf(A);
    b := evalf(B);
    while abs(b-a)>delta do
        c := (a+b)/2:
        if sign(f(c))=sign(f(a)) then
            a := c;
        else
            b := c;
        end if;
    end do;
    (a+b)/2;
end proc:

You should have done

bisection2(x -> x^2-2, 1, 3, 1e-7);

 

The 'exprseq' returned by ?whattype isn't an actual type, rather it is a name used by whattype to indicate that it was passed multiple arguments, i.e., an expression sequence.

The 'exprseq' returned by ?whattype isn't an actual type, rather it is a name used by whattype to indicate that it was passed multiple arguments, i.e., an expression sequence.

Actually, he has two expression sequences of sets.  They can be merged into a single list by first combining the expression sequences into a single list, then unpacking, with ?op, each set.

map(op, [L,G]);

or, in Maple 13

op~([L,G]);

 

Actually, he has two expression sequences of sets.  They can be merged into a single list by first combining the expression sequences into a single list, then unpacking, with ?op, each set.

map(op, [L,G]);

or, in Maple 13

op~([L,G]);

 

There are a few exceptions where a higher-level assignment to a variable used as an index in a subsequent call to add (etc.) causes grief.  The clearest is

(**) protect(i);
(**) add(i, i = 1..3);
Error, attempting to assign to `i` which is protected

In reflecting on this and the example of the escaped local k (shown previously), I assume that add (etc.) saves the previous value of the index variable, then uses it during its execution, and finally restores the previous value.  Here that fails because the protection prevents assigning to the variable.  In the escaped local case, what is happening, I assume, is that one thread has finished and restored k to its previous (unassigned value), then an active thread uses that unassigned value as the actual (incremented) value of k.  Probably it is more complicated than that, but it seems a useful narrative.

Another extension that could help with typical uses is to extend Threads:-Add, etc, to permit multiple indices. For example

Threads:-Add(i*j, i = 1..m, j = 1..n);

That wouldn't help here because of the use of the mul.
 

 

Good question.  Darin's response should be interesting.  Presumably that is one of the interesting (i.e. difficult) problems he gets to deal with. Rather than being purely automatic, could one pass in a set of variables that must be made distinct between threads?  Say

Threads:-Add(mul(i*k, k=1..10^4), i = 1..100, distinct = {k});

Seems kind of kludgy, but that could be more readily implemented.

Note that in the solution I propose (in "Explanation") it is crucial that i and j are explicitly declared as locals.  Otherwise shared global names are used and that leads to the same conflict. Would you expect that Threads:-Add and friends deal with those situations as well?

I like your idea that add, mul, and seq, which already use special evaluation rules, be modified to always use "super-local" variables.  Maybe you meant in the context of Threads, but that seems more generally useful, though it would break some code.  For example, one can currently do, at top-level

(**) x := i^2:
(**) add(x, i=1..3);
                                      14

That doesn't work in a procedure (and I don't recommend it anyway), so also wouldn't work with "super-locals":

(**) proc() local i,x; x := i^2; add(x, i=1..3); end proc();
                                        2
                                     3 i

 

That's an excellent problem that illustrates---if I understand it correctly---some of the traps of parallelizing code.

What appears to be happening is that the 'dummy' variables of the add and mul routines are being shared between different threads.  For example, in the call

Threads:-Add(add(k*i, i=1..5), k = 1..m)

all threads use the same i variable (the k is not a problem because it is handled by the Threads:-Add routine).  That causes conflicts.  I ran into that problem a few weeks ago while writing my own version of Threads:-Add (a learning project).  One way to avoid it is to call a local procedure that computes the product of sums.  That works because each call to the procedure uses its own local variables, so there is no sharing between threads. 

Here is a version that works

p2 := proc(A)
local k,l,m,n,fullL;
    m,n := op(1,A);
    fullL := [seq(1..n)];
    add((-1)^k*Threads:-Add(muladd(A,l,m), l in combinat:-choose(fullL,k)), k=1..n);
end proc:

muladd := proc(A,l,m)
local i,j;
    mul(add(A[i,j], j=l), i=1..m);
end proc:

M := Matrix(15,`+`): p2(M);
                                                 -416483579190482716042690560000

 Here's a simpler procedure that illustrates the same problem

simple := proc(n)
local j,k;
    Threads:-Add(add(k, k=1..n), j=1..n);
end proc:

simple(10^4);
                               500050342960 + k

Note the escaped k.  That won't always happen, but does occasionally. 

huh?  x=0 solves x*sin(x*y)=x for all y.  They are not equivalent.

huh?  x=0 solves x*sin(x*y)=x for all y.  They are not equivalent.

First 100 101 102 103 104 105 106 Last Page 102 of 195