Carl Love

Carl Love

28065 Reputation

25 Badges

13 years, 21 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

There is a recursive formula to count partitions which is fairly simple computationally but fairly weird mathematically (as recursive formulas go). Here it is coded in Maple:

restart:
nPart:= proc(n::integer)
option remember;
local k;
    if n<=0 then `if`(n=0,1,0)
    else
        -add(
            (-1)^k*(thisproc(n-k*(3*k-1)/2) + thisproc(n-k*(3*k+1)/2)),
            k= 1..floor((sqrt(24*n+1)+1)/6)
        )
    fi
end proc:

For large n, there are some more-complicated formulas that are faster. But the above isn't horribly slow for moderately large n.

I post this not because it's any better than combinat:-numbpart but because it's simple enough to have significant educational value.

You have 6 algebraic equations to fit 6 parameters. So it's very likely that the solution space is a 6 - 6 = 0-topological-dimension subset of C^6 (C = complex numbers), so the probability that there's an integer 6-tuple in that subset is infinitesimal. 

Note that you mistakenly used in the solve command.

Here are four non-hackish ways to do it. The first, unapply, requires that the recurrence have a closed-form solution. The others have no such requirement.

The second, makeproc, requires that the recurrence be expressible as a single algebraic formula. The remaining two do not have that requirement, and they allow formulas and procedures of arbitrary complexity.

The third, module, is a very general technique that's important to learn. They're like classes in C++. Both the local and export variables retain their values between invocations, thus avoiding the need for globals. Nearly all large pieces of Maple code are written as modules.

The fourth, procname/thisproc, is also important to learn. It uses a technique called memoization: building a table matching all input with its output. Memoization can also be made automatic by giving a procedure option remember. In this technique, the procedure's own name takes the place of the global.

restart:
Last:= unapply(rsolve({f(n)=4*f(n-1)+n, f(0)=1}, f(n)), n):
seq(Last(k), k= 1..10);

restart:
Last:= rsolve({f(n)=4*f(n-1)+n, f(0)=1}, f(n), makeproc):
seq(Last(k), k= 1..10);

restart:
Last:= module()
local _Last, ModuleApply:= (k::posint)-> (_Last:= 4*_Last+k);
export Init:= N-> (_Last:= N);
end module:
Last:-Init(1):
seq(Last(k), k= 1..10);

restart:
Last:= (k::posint)-> (procname(k):= 4*thisproc(k-1) + k):
Last(0):= 1:
seq(Last(k), k= 1..10);

 

It can be done like this:

Ab:= Matrix(M$2):
Ab[[1,-1], ..]:= A[[1,-1], ..]:
Ad:= A - Ab;

 

You have more unknowns than equations, so a purely numeric solution is not possible. I assume that you want to solve for x[0]x[1]y[0]y[1]? A solution in terms of t can be easily obtained with solve:

Sol:= solve({eq||(1..4)}, {x[0], x[1], y[0],  y[1]});

Then you can get numeric solutions for a given t by (for example),

evalf([allvalues(eval(Sol, t= 1.2))]);

The environment variable printlevel controls the nesting level to which the results of statements within do-, if-, and proc-blocks are displayed. See ?printlevel.

Yes, it has occasionally happened for me that the kernel crashes (with the boxed message) but the Save is grayed out on the File menu. In those cases, clicking the X to the right of the worksheet's name on the tabs bar will bring up a Save File dialog. (If there's only one worksheet shown on the tabs bar, that X is all the way to right edge of the window.) Either way that you save the worksheet, you'll get that kernel crash message again, but you can ignore it.

The kernel crash message says you should close Maple entirely. In my experience, this is totally false, and it's only necessary to close the worksheet(s) attached to the crashed kernel. Not closing the GUI will save you a lot of time, especially if you have a lot of open worksheets. I've done this hundreds of times without any problem.

Ordinarily, when you use evalb(...=...) to test equality, it only checks whether the objects are absolutely identical, that is, that they are in fact the exact same object with only one memory address. Many common objects in Maple, such as regular sets, are called immutable, and there's an ongoing process (happening in the background such that the user is generally unaware of it) called automatic simplification that ensures that only one copy of each such object is stored. On the other hand, for objects considered mutable, multiple copies of them that are mathematically equal may exist. And, as I said, ordinarily evalb(...=...will not say that such objects are equal. But, the object/module MultiSet overloads the = operator to bypass what ordinarily happens. You can see the overload procedure via

showstat(MultiSet:-`=`);

But, that overloading is not sophisticated enough (nor is there any existing mechanism in Maple by which it could be) to detect when you've constructed more-complicated structures which contain the MultiSets.

I'd guess that the issue with the file saving is just a case of read/write permissions. It seems unlikely to me that it has anything to do with the content of the file. To test that theory, try 

X:= 7;
save X, "file.m";
restart:
read "file.m":
X;

After 3 years, here is a solution using LPSolve for 10 deliveries with more than 1 truckload. (VV's combinatorial solution is for 8 deliveries. Brian's ILP solution is only of the 1-load TSP (i.e., travelling salesman problem = 1 load) subproblem.) The truck is required to return to the warehouse when it's finished each round of deliveries. The worksheet below specifically shows the solution for 2 loads (truck capacity K = 100). To increase the number of loads, set the value of lower. 
 

restart
:

#Vector of required deliveries. Index 1 is warehouse, so that entry is 0:
d:= <0, 67, 15, 2, 6, 12, 17, 2, 5, 13, 28>
:
#Number of destinations, including warehouse:
n:= numelems(d):
interface(rtablesize= n):
N:= {$1..n}: #index set, for convenience
C:= {$2..n}: #customer indices

#Truck capacity:
K:= 100
:

#min Truckloads required:
Vs:= ceil(add(d)/K);
V:= {$1..Vs}:

2

#Symmetric Matrix of costs (or distances) to drive between given locations:
Cost:= Matrix(
    n$2, shape= symmetric, scan= triangular,
    [0, op]~([
        [376, 246, 297, 599, 318, 132, 633, 478, 346, 585],
             [326, 240, 638, 457, 184, 746, 278, 314, 321],
                  [ 13, 619, 346,  45, 621, 323, 212, 376],
                       [645, 317, 128, 648, 379, 297, 337],
                            [545, 653, 482, 869, 742, 883],
                                 [440, 573, 340, 459, 638],
                                      [501, 365, 318, 406],
                                          [1045, 850, 1087],
                                                [243,  24],
                                                      [239]
    ])
);

Matrix(11, 11, {(1, 1) = 0, (1, 2) = 376, (1, 3) = 246, (1, 4) = 297, (1, 5) = 599, (1, 6) = 318, (1, 7) = 132, (1, 8) = 633, (1, 9) = 478, (1, 10) = 346, (1, 11) = 585, (2, 1) = 376, (2, 2) = 0, (2, 3) = 326, (2, 4) = 240, (2, 5) = 638, (2, 6) = 457, (2, 7) = 184, (2, 8) = 746, (2, 9) = 278, (2, 10) = 314, (2, 11) = 321, (3, 1) = 246, (3, 2) = 326, (3, 3) = 0, (3, 4) = 13, (3, 5) = 619, (3, 6) = 346, (3, 7) = 45, (3, 8) = 621, (3, 9) = 323, (3, 10) = 212, (3, 11) = 376, (4, 1) = 297, (4, 2) = 240, (4, 3) = 13, (4, 4) = 0, (4, 5) = 645, (4, 6) = 317, (4, 7) = 128, (4, 8) = 648, (4, 9) = 379, (4, 10) = 297, (4, 11) = 337, (5, 1) = 599, (5, 2) = 638, (5, 3) = 619, (5, 4) = 645, (5, 5) = 0, (5, 6) = 545, (5, 7) = 653, (5, 8) = 482, (5, 9) = 869, (5, 10) = 742, (5, 11) = 883, (6, 1) = 318, (6, 2) = 457, (6, 3) = 346, (6, 4) = 317, (6, 5) = 545, (6, 6) = 0, (6, 7) = 440, (6, 8) = 573, (6, 9) = 340, (6, 10) = 459, (6, 11) = 638, (7, 1) = 132, (7, 2) = 184, (7, 3) = 45, (7, 4) = 128, (7, 5) = 653, (7, 6) = 440, (7, 7) = 0, (7, 8) = 501, (7, 9) = 365, (7, 10) = 318, (7, 11) = 406, (8, 1) = 633, (8, 2) = 746, (8, 3) = 621, (8, 4) = 648, (8, 5) = 482, (8, 6) = 573, (8, 7) = 501, (8, 8) = 0, (8, 9) = 1045, (8, 10) = 850, (8, 11) = 1087, (9, 1) = 478, (9, 2) = 278, (9, 3) = 323, (9, 4) = 379, (9, 5) = 869, (9, 6) = 340, (9, 7) = 365, (9, 8) = 1045, (9, 9) = 0, (9, 10) = 243, (9, 11) = 24, (10, 1) = 346, (10, 2) = 314, (10, 3) = 212, (10, 4) = 297, (10, 5) = 742, (10, 6) = 459, (10, 7) = 318, (10, 8) = 850, (10, 9) = 243, (10, 10) = 0, (10, 11) = 239, (11, 1) = 585, (11, 2) = 321, (11, 3) = 376, (11, 4) = 337, (11, 5) = 883, (11, 6) = 638, (11, 7) = 406, (11, 8) = 1087, (11, 9) = 24, (11, 10) = 239, (11, 11) = 0})

#Check and ENFORCE triangle inequality:
for T in combinat:-choose(n,3) do
    if (S:= Cost[T[1],T[2]] + Cost[T[2],T[3]]) < Cost[T[1],T[3]] then
        printf("Triangle inequality violation at %a corrected.\n", T);
        Cost[T[1],T[3]]:= S #Enforcement
    fi
od:  

Triangle inequality violation at [1, 3, 4] corrected.
Triangle inequality violation at [1, 7, 11] corrected.
Triangle inequality violation at [1, 9, 11] corrected.
Triangle inequality violation at [2, 7, 8] corrected.

Triangle inequality violation at [2, 9, 11] corrected.
Triangle inequality violation at [3, 4, 6] corrected.
Triangle inequality violation at [3, 4, 11] corrected.
Triangle inequality violation at [3, 7, 8] corrected.
Triangle inequality violation at [3, 9, 11] corrected.
Triangle inequality violation at [4, 7, 8] corrected.
Triangle inequality violation at [6, 9, 11] corrected.
Triangle inequality violation at [7, 9, 11] corrected.
Triangle inequality violation at [8, 9, 11] corrected.

#View corrected matrix:
Cost;

Matrix(11, 11, {(1, 1) = 0, (1, 2) = 376, (1, 3) = 246, (1, 4) = 259, (1, 5) = 599, (1, 6) = 318, (1, 7) = 132, (1, 8) = 633, (1, 9) = 478, (1, 10) = 346, (1, 11) = 502, (2, 1) = 376, (2, 2) = 0, (2, 3) = 326, (2, 4) = 240, (2, 5) = 638, (2, 6) = 457, (2, 7) = 184, (2, 8) = 685, (2, 9) = 278, (2, 10) = 314, (2, 11) = 302, (3, 1) = 246, (3, 2) = 326, (3, 3) = 0, (3, 4) = 13, (3, 5) = 619, (3, 6) = 330, (3, 7) = 45, (3, 8) = 546, (3, 9) = 323, (3, 10) = 212, (3, 11) = 347, (4, 1) = 259, (4, 2) = 240, (4, 3) = 13, (4, 4) = 0, (4, 5) = 645, (4, 6) = 317, (4, 7) = 128, (4, 8) = 629, (4, 9) = 379, (4, 10) = 297, (4, 11) = 337, (5, 1) = 599, (5, 2) = 638, (5, 3) = 619, (5, 4) = 645, (5, 5) = 0, (5, 6) = 545, (5, 7) = 653, (5, 8) = 482, (5, 9) = 869, (5, 10) = 742, (5, 11) = 883, (6, 1) = 318, (6, 2) = 457, (6, 3) = 330, (6, 4) = 317, (6, 5) = 545, (6, 6) = 0, (6, 7) = 440, (6, 8) = 573, (6, 9) = 340, (6, 10) = 459, (6, 11) = 364, (7, 1) = 132, (7, 2) = 184, (7, 3) = 45, (7, 4) = 128, (7, 5) = 653, (7, 6) = 440, (7, 7) = 0, (7, 8) = 501, (7, 9) = 365, (7, 10) = 318, (7, 11) = 389, (8, 1) = 633, (8, 2) = 685, (8, 3) = 546, (8, 4) = 629, (8, 5) = 482, (8, 6) = 573, (8, 7) = 501, (8, 8) = 0, (8, 9) = 1045, (8, 10) = 850, (8, 11) = 1069, (9, 1) = 478, (9, 2) = 278, (9, 3) = 323, (9, 4) = 379, (9, 5) = 869, (9, 6) = 340, (9, 7) = 365, (9, 8) = 1045, (9, 9) = 0, (9, 10) = 243, (9, 11) = 24, (10, 1) = 346, (10, 2) = 314, (10, 3) = 212, (10, 4) = 297, (10, 5) = 742, (10, 6) = 459, (10, 7) = 318, (10, 8) = 850, (10, 9) = 243, (10, 10) = 0, (10, 11) = 239, (11, 1) = 502, (11, 2) = 302, (11, 3) = 347, (11, 4) = 337, (11, 5) = 883, (11, 6) = 364, (11, 7) = 389, (11, 8) = 1069, (11, 9) = 24, (11, 10) = 239, (11, 11) = 0})

#
# We define the decision variables thus:
#    x[i,j,v] = 1, if truckload v goes from destination i to j;
#             = 0, otherwise.
#

#Objective:   (Eq 2.2)
#   minimize:
Obj:= add(add(add(Cost[i,j]*x[i,j,v], i= N), j= N), v= V)
:

###############
# Constraints #
###############
Con:= table()
:
# Each truckload does not exceed the capacity: (Eq 2.4)
#
Con[capacity]:=
    seq(add(d[j]*add(x[i,j,v], i= N minus {j}), j= N) <= K, v= V)
:

# Each delivery is made exactly once:  (Eq 2.3)
#
Con[delivery]:=
    seq(add(add(x[i,j,v], j= N minus {i}), v= V) = 1, i= C)
:

# Trucks start at warehouse:  (Eq 2.5)
#
Con[start]:=
    seq(add(x[1,j,v], j= C) = 1, v= V)
:

# Trucks return to warehouse: (Carl's idea)
#
Con[`return`]:=
    seq(add(x[i,1,v], i= C) = 1, v= V)
:

# Continuity: (Eq 2.6)
#
Con[continuity]:=
    seq(
        seq(
            add(x[i,k,v], i= N minus {k})
                = add(x[k,j,v], j= N minus {k}),
            k= C
        ),
        v= V
    )
:

# Avoid subtours:  (VV's idea) (u's are artificial variables.)
#
Con[antisubtour]:=
    seq(
        seq(
            seq(
                u[i,v]-u[j,v]+n*x[i,j,v] <= n-1,
                i= N minus {1,j}
            ),
            j= C
        ),
        v= V
    )
:

BV:= binaryvariables
    = {seq(seq(seq(x[i,j,v], i= N minus {j}), j= N), v= V)}
:

infolevel[Optimization]:= 5:
Sol:= CodeTools:-Usage(
    Optimization:-LPSolve(
        Obj, [entries(Con, nolist)],
        BV, assume= nonnegative, feasibilitytolerance= 1e-5
    )
):

LPSolve: calling ILP solver

LPSolve: number of problem variables 240
LPSolve: number of general linear constraints 216
LPSolve: feasibility tolerance set to 0.1e-4
LPSolve: integer tolerance set to 0.10e-4
LPSolve: infinite bound set to 0.10e21
LPSolve: iteration limit set to 2280
LPSolve: depth limit set to 360
LPSolve: node limit set to 0

memory used=5.86MiB, alloc change=2.76MiB, cpu time=29.52s, real time=29.68s, gc time=0ns

Stops:= {lhs~(select(e-> rhs(e)::1, Sol[2]))[]}:

Tour:= table():
for X in Stops do
    Tour[op(3,X)][op(1,X)]:= op(2,X)
od:

for v in V do
    A:= Array(1..1, [1]);
    for i do A,= (r:=  Tour[v][A[i]]) until r=1;
    Tour[v]:= seq(A)
od:

Routes:= [entries(Tour)];

[[1, 2, 10, 1], [1, 6, 9, 11, 4, 3, 7, 8, 5, 1]]

COST = Sol[1];

COST = 3695

 


 

Download VehicleRoutingILP.mw

Do you mean something like this?

unapply([for k to 3 do x^k od], x);

This uses syntax new to Maple 2019. The example above would be better done with seq, but there are more-elaborate examples where a do-loop (or for-loop) works and seq doesn't. There are also cases where seq is the less-efficient option.

I would've given the same answer as Joe and Acer if I'd gotten to it in time. I'll add that the same level of type checking can also be applied to the return value of a procedure:

P:= proc()::posint; 0 end proc:  

P(); #error when assertlevel is set to 2

But mostly I want to draw your attention to the distinction between types and properties because I think that you may be referring to the latter. Types are inherent to data structures (i.e., anything that you can assign to a variable): Given any structure and any type ttype(x,t) should return true or false; nebulous or ambiguous cases should be avoided when coding t.

Properties, on the other hand, can be given to symbolic variables (i.e., names) with assume and related commands. Some "property-aware" commands will respect the properties, but most commands will ignore them.

Maple's type system is highly developed, but the property system is nascent and has been so AFAIK since created by Gaston Gonnet, et al, many years ago. My guess is that making significant improvements to the system would require redesign of Maple from the ground up. I believe that the strongly typed "free" CAS named Axiom has something like this implemented. Maybe it'd be possible to get Maple and Axiom to interact?

Use:

cat("", V, ".xls")

The plot3d command needs more quotes. This works:

plot3d(['f(x,y,'g(x, y)')'], x=0..1,y=0..1)

Use

f~([$0..5])

If you put immediately before an operator, it makes it inert. So, you can do a %/ b * c, and it'll display like you want. The inertness will remain until you use the value command.

Unfortunately, you do need to use some special syntax to get things to display as you want.

If an appears in the first column of input, it's a shell escape, equivalent to a system command. See ?system.

First 118 119 120 121 122 123 124 Last Page 120 of 395