Carl Love

Carl Love

28015 Reputation

25 Badges

12 years, 298 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Erik,

Here's a complete example of using the recursive procedure and generating the output with evalhf. I used the sum-of-squares for the first real-valued evaluation function, and the sum-of-cubes for the second.

Combi:= module()
export
    boxContent,  #A vector whose i'th entry holds the number of marbles in the i'th box
    Output       #Matrix of Output value
;
local
    Combi:= proc(
        k::integer, Combi:= module()
        boxContent::Vector,
        activeBox::integer,
        marblesRemaining::integer,
        Output::Matrix,
        Index::Vector
    )       
    local i::integer;
        
        if activeBox < k then
            for i from 0 to marblesRemaining do
                boxContent[activeBox]:= i;
                thisproc(k, boxContent, activeBox+1, marblesRemaining-i, Output, Index)
            end do
        else
            boxContent[activeBox]:= marblesRemaining;
            Index[1]:= Index[1] + 1;
            #Something:
            Output[1, Index[1]]:= add(boxContent[i]^2, i= 1..k);
            Output[2, Index[1]]:= add(boxContent[i]^3, i= 1..k);        
        end if;
    end proc,

    ModuleApply:= proc(n::posint, k::posint)
    local
        Index,
        boxContent  #A vector whose i'th entry holds the number of marbles in the i'th box
    ;
        boxContent:= Vector(k, datatype= float[8]);
        Output:= Matrix(2, binomial(n+k-1, k-1), datatype= float[8]);
        Index:= Vector(1);
        evalhf(Combi(k, boxContent, 1, n, Output, Index));
        [][]
    end proc
;
end module;

And the worksheet:


A recursive procedure

restart

with(Statistics):

n := 40:

k := 6:

 

module () local Combi, ModuleApply; export boxContent, Output; end module

NULL

B := CodeTools[Usage](Combi(n, k));

memory used=18.65MiB, alloc change=18.64MiB, cpu time=3.45s, real time=3.45s

Combi:-Output[1, 1221759];

HFloat(1600.0)

Combi:-Output[2, 1221759];

HFloat(64000.0)

``


Download A_recursive_procedu.mw

From your worksheet, I think that you've already noticed that this is a bug related to the strange way that you created V, which is equivalent to this shortened form:

restart:
v1:= <a>: v2:= <b>:
V:= Vector(2, [v1,v2]):
subs(a=5, V);

The workaround is to create V like this:

V:= <v1,v2>:
subs(a=5, V);

RootPlot:= (p::polynom)->
     plot(map([Re,Im], [fsolve](p, indets(p,name)[], complex)), style= point,  _rest):

RootPlot(1+randpoly(x, degree= 2^7), color= green, symbol= diamond);

The command StringTools:-ParseTime (see ?StringTools,ParseTime ) can be used on the timestamps, like this (ignoring the time zone):

restart:
timestamp:= "2013-08-13T00:29:24+0000": #Note the quotes
timerecord:= StringTools:-ParseTime("%Y-%m-%dT%T", timestamp):

This produces a structure called a ?Record which holds then numeric values for year, month, etc.

baseyear:= 2012:  basemonth:= 9:  #For example
CountMonths:= TR-> 12*(TR:-year - baseyear) + TR:-month - basemonth:
CountMonths(timerecord);
                               11
Now this gives you a way to measure distance on the horizontal axis.

The command ListTools:-Classify can be used to count the photos, texts, etc., for each month (see ?ListTools,Classify ).

See ?plot,tickmarks for information about changing the tickmarks on the horizontal axis.

See ?plots,polygonplot to plot the colored rectangles.

 

The problem is that you did not declare an Array, so gg became a table. But, gg cannot be a Vector (or a list) because it is indexed from 0 instead of 1. It has to be an Array or a table.

You can create an Array and fill it in one step:

gg:= Array(0..7, i-> t-> t^i):

If you want to apply all the functions in gg at once, you can do this (I am using 2 as the argument):

gg2:= map(`?()`, gg, [2]);

gg2[5];
                          32

S:= {[-2, 1, 1, -2, 1, 1], [-2, 1, 1, 2, -1, -1], [0, -1, 1, 0, 1, -1], [0, 1, -1, 0, -1, 1]}:

Vectors u and v are multiples iff (u.v)*u = (u.u)*v. Proof of this is left as an exercise.

IsMultiple:= (u::~Vector, v::~Vector)-> ArrayTools:-IsZero(u.v*u - u.u*v):

If we exclude the zero vector, which is a multiple of everything, IsMultiple is an equivalence relation on the remaining vectors. ListTools:-Categorize will partition a list into equivalence classes.

S1:= S minus {[0$nops(S[1])]};   
L:= ListTools:-Categorize(IsMultiple, [S1[]]);
       [[-2, 1, 1, -2, 1, 1]], [[-2, 1, 1, 2, -1, -1]],
         [[0, -1, 1, 0, 1, -1], [0, 1, -1, 0, -1, 1]]

Select a representative (the first element) of each equivalence class.

map(`?[]`, {L}, [1]);
{[-2, 1, 1, -2, 1, 1], [-2, 1, 1, 2, -1, -1], [0, -1, 1, 0, 1, -1]}

The answer is, briefly, no, Maple cannot work with matrices and vectors of arbitrary, symbolic dimension. (Considering the number of questions that I see about it here, this could be one of the next frontiers for symbolic computation.) If you can deduce (with mindpower) a recurrence along a single dimension, you might be able to solve it with rsolve. Multi-dimensional recurrences are another frontier.

This is of no consequence in your program currently, but the statement

r:= Vector[(N1+1)*N2+1];

in your procedure main doesn't do anything useful. If you did want to pre-allocate a Vector, the statement should be with ( ) instead of [ ].

subsindets(
     test1r2,
     `^`,
     x-> `if`(op(2,x) = 1/2,`if`(op(1,x)::`+`, `Math.Sqrt`(op(1,x)), x), `Math.Pow`(op(x)))
);

Answering your second question, there is the command interface(prompt= ...). For example, I have the following line in my initialization file:

interface(prompt= "(**)");

Erik wrote:

I am working on a problem, which basicly is dealing with putting n identical marbles in k different boxes. I need to run through all possible combinations and do something with it.

Your "combinations" are a combinatorial structure called "compositions".

I have heard that a recursive procedure sometimes can be transformed into an iterative procedure, which will maybe make it execute faster.

That is for a small but important subset of recursions called "tail recursion". Your procedure is not that type. An example of tail recursion is

Factorial:= (n::nonnegint)-> `if`(n=0, 1, n*Factorial(n-1));

The issue with recursion is usually memory more than time.

I can't find a way to do it iteratively....

I found two pre-coded iterative ways in Maple. The first uses a very old package called combstruct. The second uses Joe Riel's (excellent) Iterator package available for free download at the Maple Applications Center here. The problem is discussed in Donald Knuth's (excellent and long-awaited) The Art of Computer Programming: Volume 4. But these methods are much slower than your recursive method. The reason to use an iterative method is that with a recursive method you essentially have to hold all combinations in memory at the same time. This is no problem for the binomial(45,5) ~ 1.3 M combinations in your example. But if you just increase your number of boxes to 10, then you have binomial(49,9) ~ 2 G, each one being an Array with 10 elements, and it's not possible in the memory of most computers.

With Joe's Iterator package, the code for your procedure is simply:

Combi:= proc(n::posint, k::posint)
local boxContent;
     for boxContent in Iterator:-BoundedComposition([n$k], n) do
          # Do something here
     end do
end proc;

With combstruct, it is

Combi:= proc(n::posint, k::posint)
local
     Iter:= combstruct[iterstructs](Composition(n+k), size= k),
     NV:= Iter[nextvalue],
     boxContent
;
     while not Iter[finished] do
          boxContent:= NV() -~ 1;
          # Do something here
     end do
end proc;

If we need to stay with this recursive procedure, then is there a way to optimize the code?

Like I said, your procedure is already quite efficient. I can make a factor-of-five improvement by using evalhf, but this will only be possible if your "# Do something here" is doable within evalhf, which is a very restricted environment. Here is your procedure done with evalhf:

Combi:= module()
export boxContent;  #A vector whose i'th entry holds the number of marbles in the i'th box

local
    Combi:= proc(k::float, boxContent::Vector, activeBox::float, marblesRemaining::float)       
    local i::integer;
        
        if activeBox < k then
            for i from 0 to marblesRemaining do
                boxContent[activeBox]:= i;
                thisproc(k, boxContent, activeBox+1, marblesRemaining-i)
            end do
        else
            boxContent[activeBox]:= marblesRemaining;
            #Do something here        
        end if;
    end proc,

    ModuleApply:= proc(n::posint, k::posint)
        boxContent:= Vector(k, datatype= float[8]);
        evalhf(Combi(k,boxContent,1,n))
    end proc
;
end module;

The call for all three procedures is simply

Combi(n, k)

I was able to do it in under a second by using simplify with side relations (see ?simplify,siderels ). I refer to it here by its syntax skeleton simplify(..., [...]). Unlike regular simplify---which may use heuristic ideas of what is simpler---simplify(..., [...]) uses a pure polynomial algebraic algorithmic approach. 

Replacing a common subexpression with a new variable is a type of side relation. Here are the substitutions that you wanted to apply:

There are a number of factors that can make a vast difference in the amount of time and memory that the algorithm uses. These factors are referred in the help page as "monomial orderings", but the help page fails to mention just how important these are. The naive approach to your problem would be

simplify(kappa, [newpar]);

I let this run several minutes until I was nearly out of memory (it was using about 4Gig) before I killed it. Then I tried applying the side relations individually from most complicated to least (this last point is very important), like this:

newkappa:= kappa:
for k from nops([newpar]) by -1 to 1 do
     newkappa:= simplify(newkappa, [newpar[k]])
end do:
newkappa;

This works nearly instantaneously (0.140 secs) and achieves a great simplification: The only variables in the result are the s's. I also did the C=4, K=7 case in 0.656 secs. So, there you go: I took it from several days of computation to under a second.

Here's the modified worksheet:

sC.mw

In the line before the else, there appears #Iota.... The # comments out the rest of the line, so it unbalances the parentheses.

Your erroneous result seems to be related to a bad simplification of factorials. After evaluating at x=1 and r=1, your summand has the subexpression

ex:= (1+k)!/k!^2/(2-k)!:

which simplifies as

simplify(ex);

Obviously this is wrong for any integer k. A workaround is assume(k::nonnegint).

simplify(ex) assuming k::nonnegint;

To get rid of all of the RootOfs at once, you could do

evalindets(newpar, RootOf, x-> allvalues(x)[1]);

First method (my preference):  second:= eval(t, test[2]);

Second method: second:= rhs(test[2][]) or rhs(test[2][1]) or op([2,1,2], test)

Unfortunately, there is no analog to the first method for selecting the first value; it only works on expressions of the form (variable = value).

First 355 356 357 358 359 360 361 Last Page 357 of 394