Joe Riel

9530 Reputation

23 Badges

20 years, 28 days

MaplePrimes Activity


These are answers submitted by Joe Riel

Try ?PolynomialTools[CoefficientVector].

Here's a module that you can use to analyze the Parrondo paradox.

Parrondo := module()
option package;
export ExpectedValue
    ,  Pstationary
    ,  ProbOfEvents
    ,  A, B;
    ;
global e;

    # Define the state-transition matrices.
    # Pij is the probability of going from
    # state i to state j.  
    
    A := Matrix(3, [[   0   , 1/2-e , 1/2+e ],
                    [ 1/2+e ,   0   , 1/2-e ],
                    [ 1/2-e , 1/2+e ,   0   ]]
               );
    
    B := Matrix(3, [[   0   , 1/10-e , 9/10+e ],
                    [ 1/4+e ,   0    , 3/4-e  ],
                    [ 3/4-e , 1/4+e  ,   0    ]]
               );
    
    Pstationary := proc( P :: list(Matrix) )
    local lambda,M,V,y,Y,pos;
    description "Compute the stationary probability vector";
        # M is the transition matrix of the system.
        M := foldl(`.`, P[]);
        # The stationary probability vector is a row vector 
        # that satisfies Ps.M = Ps.  Use Eigenvectors
        # with the transpose of M to compute the tranpose of Ps.
        (lambda,V) := LinearAlgebra:-Eigenvectors(M^%T);
        if not member(1, convert(lambda,list), 'pos') then
            error "cannot find unity eigenvalue";
        end if;
        Y := V[..,pos];
        # Normalize the vector.
        return normal~(Y/add(y, y in Y));
    end proc;
    
    ExpectedValue := proc( P :: list(Matrix) )
    local i,k,bits,evts,val,E,n,Ps;
    description "Compute the expected value";
        Ps := Pstationary(P);
        n := nops(P);
        E := 0;
        for k from 0 to 2^n-1 do
            # Sum value of each possible sequence of evts.
            bits := Bits:-Split(k,':-bits'=n);
            evts := subs(0=-1, bits);
            val := add(e, e in evts);
            if val <> 0 then
                E := E + val*add(Ps[i]*ProbOfEvents(i, evts, P), i = 1..3)
            end if;
        end do;
        return normal(E);
    end proc;
    
    ProbOfEvents := proc(i0::{1,2,3}, evts::list({-1,+1}), P :: list(Matrix))
    description "Compute the probability of a given sequence of events, starting at state i0";
    local i,j,k,p;
        i := i0;
        p := 1;
        for k to nops(evts) do
            j := i + evts[k];
            j := 1 + modp(j-1, 3);
            p := p*P[k][i,j];
            i := j;
        end do;
        return p;
    end proc;
    
end module:

For example, to compute the expected (stationary) value of the game A, A, B, B, do

with(Parrondo):
E := ExpectedValue([A,A,B,B]);
eval(E, e=0);
                                           16/163
plot(E, e=0..0.02);

Be careful with has(M,1).  It is not documented to work only on the entries of rtables, though it appears to.  That is, from the help page one might conclude that

has(Matrix(1,1,[[3]]) 

should return true because one of its subexpressions from ?op does:

op(Matrix(1,[3]));        
                1, 1, {(1, 1) = 3}, datatype = anything, storage = rectangular, order = Fortran_order, shape = []

That is not the case, but that doesn't necessarily mean it does what you want.  As an example (possibly not relevant):

has(1/x, -1);
                                   true

 

Earlier this week I wrote sorting with a permutation list, which shows how to do just that. That is, it generates a permutation that can then be applied to x to get y. Applying it here

x := [1,5,3,7]:
P := map(attributes, sort([seq(setattribute(evalf(x[i]),i), i = 1..nops(x))]));
                         P := [1, 3, 2, 4]
[seq(x[P[i]], i = 1..4)];
                         [1, 3, 5, 7]

You could avoid the issue by using rationals, say 1/10 or 1/20.

If Maple would display an overloaded operator in infix rather than prefix form, then it would be easy:

delayadd := module()
option package;
export `+`;
end module:

with(delayadd):
1+2+3+4;
                           +(+(+(1, 2), 3), 4)

 

You could do

delayadd := module()
option package;
export `+`;
    `+` := proc() :-`+`(map(``, [_passed])[]) end proc;
end module:

with(delayadd):
1+2+3+4;
                         ( ( (1) +  (2)) +  (3)) +  (4)

but the resulting nesting isn't very nice

You can do substantially better, but with a bit of work.  First, let's compute an expression for the number of valid lists.  Given the integers 1..n, we know that the minimum allowable k (number of zeros) is n-1.  Let z be the number of excess zeros, that is, z = k-n+1.  For the integers 1..n-1, each must be followed by a zero.  So think of the integers as the pairs (1,0), (2,0), ... (n-1,0) and the singleton n.  We have z zeros remaining to place among these n objects.  This is equivalent to choosing z objects from n+z objects, or binomial(n+z,z) = binomial(k+1, k-n+1).  So for the n=6, k=8 example, the number of lists is binomial(9,3) = 84.

To make this sort of thing really fast, in Maple, consider compiling the code.  Furthermore, to reduce the memory usage and gain speed, generate and operate on each output list one at at a time, rather than generating all the lists at once.  More accurately, rather than generating a list of lists, replace, in a loop, the contents of a mutable data structure (an rtable).

I have written a module that I hope to release soon (still needs some cleanup) that creates fast (compilable) iterators for common finite structures. I'll demonstrate it here. The procedure that generates each of the structures you want is

myp1 := proc(n,k)
local z, P,A,Z,W;
    z := k-n+1; # excess zeros
    (P,A,Z) := Iterator:-Combination(n+z,z,'compile');
    W := Array(1..n+k, 'datatype' = integer[4]);  # this Array stores the desired result
    while P(A) do  
        LocToComb1(n,z,Z,W);
         printf("%d\n", W);  # normally one would do something useful with the generated Array at this point
    end do;
end proc:

Procedure to convert the Array Z, which contains the logical positions of the excess zeros, to the final form, which is stored in W.

LocToComb1 := proc(n,z, Z :: Array(datatype=integer[4]), W :: Array(datatype=integer[4]))
option autocompile;
local i,j,k,p,posz;
    p := 1;  # index into W
    j := 1;  # index into Z
    k := 1;  # integer index
    posz := Z[j]+1;
    for i to n+z do
        if i = posz then
            # insert excess zero
            W[p] := 0;
            p := p+1;
            if j < z then
                j := j+1;
                posz := Z[j]+1;
            else
                # no more excess zeros; finish up integers
                for k from k to n-1 do
                    W[p] := k;
                    W[p+1] := 0;
                    p := p+2;
                end do;
                W[p] := n;
                break;
            end if;
        else
            # insert integer
            W[p] := k;
            if k < n then
                # inset paired zero
                W[p+1] := 0;
                p := p+2;
                k := k+1;
            else
                # no more integers; finish up
                for p from p+1 to 2*n-1+z do
                    W[p] := 0
                end do;
                break;
            end if;
        end if;
    end do;
    return 0; # not used other than to define return type for compiler
end proc:

Here is an example,
 

myp1(3,4);
0 0 1 0 2 0 3
0 1 0 0 2 0 3
1 0 0 0 2 0 3
0 1 0 2 0 0 3
1 0 0 2 0 0 3
1 0 2 0 0 0 3
0 1 0 2 0 3 0
1 0 0 2 0 3 0
1 0 2 0 0 3 0
1 0 2 0 3 0 0

How fast is it?

time(myp1(6,18)); # with printing disabled
                                       0.140

Contrast this with your procedure

time(myp2(6,18));
                                      5.256

Followup

A simpler way to do this is to generate all permutations of n ones and z zeros, then substitute the k'th 1 with the integer pair (k,0), omitting the zero if at the last position.  This can be done with

PermToPattern := proc(n :: integer, z :: integer
                      , Z :: Array(datatype=integer[4])
                      , W :: Array(datatype=integer[4])
                     )
local i,j,p;
option autocompile;
    j := 1;
    p := 1;
    for i to n+z do
        if Z[i] = 0 then
            W[p] := 0;
            p := p+1
        else
            W[p] := j;
            j := j+1;
            if j > n then
                # done with integers; fill remaining 0's
                for p from p+1 to 2*n+z-1 do
                    W[p] := 0;
                end do;
                return 0;
            end if;
            W[p+1] := 0;
            p := p+2;
        end if;
    end do;
    return 0;
end proc;


myp4 := proc(n,k)
local z, P,A,Z,W;
    z := k-n+1;
    (P,A,Z) := Iterator:-Permute([1$n,0$z],'compile');
    W := Array(1..n+k, 'datatype' = integer[4]);
    while P(A) do
        PermToPattern(n,z,Z,W);
        printf("%d\n", W);
    end do;
end proc:

This is not any faster to run, but it was easier to write the procedure that transforms each permutation into the desired pattern.

time(myp4(6,18));
                                      0.168

A somewhat more general approach, say if B included sin(u) into which you wanted to subsitute, is

frontend(curry(subs, A), [B], [{Not(specfunc(anything, {f,D(f)}))},{}]);

 

The solution that Doug gave, in the original post, was close to optimal Maple technique; it was short, efficient, and clear, once you understand how lists can be indexed from either end (which he explained).  But, if you don't understand it, or cannot explain it adequately, then you should try something else. 

Here's another approach. I won't write any code, so you cannot be accused of cutting and pasting.  Split the sorted list into two equal length lists (if there is an odd number of elements, the second list should be one smaller than the first. You can use ?op or list indexing to do so.  Use ?ListTools[Reverse] to reverse the order of the second list.  Then use ?ListTools[Interleave] to interleave the contents of the two lists.  That should give you what you want.  It isn't as nice as Doug's solution, but if you are able to implement it, then it should suffice.  As an example, you would do

[1, 2, 3, 4, 5, 6, 7]     # start with this
[1, 2, 3, 4], [5, 6, 7]   # split
[1, 2, 3, 4], [7, 6, 5]   # reverse second list
[1, 7, 2, 6, 3, 5, 4]     # interleave results

 

Here's one approach.  It runs in about 0.1 seconds and is reasonably short.  Each line of output consists of the rows in which to place a queen.

queens := module()

export ModuleApply;

local guess
    , init
    , place
    , Solve
    , N,cnt,board
    ;

    ModuleApply := Solve;

    init := proc(n)
        N := n;
        board := Array(1..N,1..N);
        cnt := 0;
        NULL;
    end proc;

    Solve := proc(n)
        init(n);
        guess(1);
        return cnt;
    end proc;

    guess := proc(col)
    local row;
        for row to N do
            if board[row,col] = 0 then
                place(row,col,+1);
                if col = N then
                    cnt := cnt+1;
                    print(_rest,row);
                else
                    procname(col+1,_rest,row)
                end if;
                place(row,col,-1);
            end if;
        end do;
    end proc;

    place := proc(row,col,val)
    local i;
        board[row,col] := -val-1;
        for i from col+1 to N do
            board[row,i] := board[row,i]+val;
        end do;
        for i to min(row-1,N-col) do
            board[row-i,col+i] := board[row-i,col+i] + val;
        end do;
        for i to min(N-row,N-col) do
            board[row+i,col+i] := board[row+i,col+i] + val;
        end do;
    end proc;

end module:
queens(8);

 

You might upload the file so we could take a look at it. If large, consider zipping it first.

I don't understand the difficulty.  You can put as many commands as you want in the loop.

Here's what I meant to post,

with(LinearAlgebra):
M := <<m+lambda|m>,<m|2*m>>:
ans := solve(Determinant(M),[m]);
                       ans := [[m = 0], [m = -2 lambda]]


M2 := map(subs, ans, M);
                       [lambda    0]  [ -lambda     -2 lambda]
                M2 := [[           ], [                      ]]
                       [  0       0]  [-2 lambda    -4 lambda]



for l from 1 to 2 by 0.1 do
    lambda=l, subs(lambda = l, M2);
end do;
                                   [1    0]  [-1    -2]
                      lambda = 1, [[      ], [        ]]
                                   [0    0]  [-2    -4]

                                 [1.1    0]  [-1.1    -2.2]
                  lambda = 1.1, [[        ], [            ]]
                                 [ 0     0]  [-2.2    -4.4]

...

Use ?evalf to evaluate them to floating-points.  You also might consider using ?fsolve rather than solve.

There are a few problems with your computation.  First, you solve for m, using the determinant of M, then plug into MM, which contains m[i] rather than m; that doesn't make sense to me.  Second, at the exit of the for loop, lambda is assigned a numerical value (it is the loop counter).  I'm thinking you don't want that.  Here is what I might do,
 

with(LinearAlgebra):
M := <<m+1|m>,<m|2*m>>:
ans := solve(Determinant(M),[m]);
                          ans := [[m = 0], [m = -2]]


for ans in ans do
    subs(ans, M);
end do;
                                   [1    0]
                                   [      ]
                                   [0    0]

                                  [-1    -2]
                                  [        ]
                                  [-2    -4]

Note the bit of trickery using ans as the loop counter. Normally I wouldn't do that, but here it illustrates that it can be done, the original value is independent of the loop value, though at the conclusion ans is modified (equal to the last value, [m=-2]). Instead of a loop, you could do

map(subs, ans, M);
                             [1    0]  [-1    -2]
                            [[      ], [        ]]
                             [0    0]  [-2    -4]

First 73 74 75 76 77 78 79 Last Page 75 of 114