Robert's solution is excellent and probably exactly what the asker wanted. However, I can't help but suggest a version that's substantially more efficient for very big matrices (and, unavoidably, a few lines longer):

listshuffle := (mat :: Matrix) -> Matrix(op(1, mat), combinat['randperm'](convert(mat, 'list')));
statsshuffle := proc(mat :: Matrix)
local v :: Vector;
  v := ArrayTools:-Alias(mat, [`*`(op(1, mat))]);
  v := Statistics:-Shuffle(v);
  return ArrayTools:-Alias(v, [op(1, mat)]);
end proc;

test := Matrix(1000, 1000, symbol=a);
CodeTools:-Usage(listshuffle(test));
CodeTools:-Usage(statsshuffle(test));

On my machine, the listshuffle takes 2.89 seconds and 235 MB memory, whereas the statsshuffle takes 0.11 seconds and only 7.7 MB. This is because in statsshuffle:
1) there is only one occasion where the data is copied, and no conversions (ArrayTools:-Alias does not actually convert the data in the rtable that you pass to it, it just puts a different face on it),
2) that one copy occurs in Statistics:-Shuffle, which does its work in fast external code, and works with rtables only.
listshuffle, on the other hand, first converts the big matrix to a list (memory-intensive), then computes a random permutation of the integers in external code (fast), and then permutes the list according to that random permutation (slow, involving a lot of random memory access). It then recreates the Matrix (memory-intensive again).

To wrap this up, if you need to do something like this in a program, consider using tools that do not require you to convert your data type - and see if you can subvert ArrayTools:-Alias to your purposes!

Hope this helps,

Erik Postma
Maplesoft.

PS: I just redid the test with a datatype=float[8] Matrix, and while the timings and memory usage are no different, statsshuffle has the advantage that the Matrix returned automatically has datatype=float[8]. For listshuffle we would need to explicitly tell Maple to reuse the same rtable attributes.

I can't believe I didn't check wikipedia: http://en.wikipedia.org/wiki/Higher-order_singular_value_decomposition explains that there are two notions of SVD for higher-order tensors, and the one I need is called candecomp-PARAFAC. I guess I'll be implementing that some time next week...

Erik.

What a short bike ride can't do for you (from work to home). Here is at least an approximation algorithm, valid for general kth order tensors v.

Let u be the best rank-one approximation to v. If v = u, stop, otherwise approximate v - u recursively.

Now the problem has been reduced to finding the best rank one approximation to a general kth order tensor v. I think we can do this as follows:

Add the entries along one dimension of the tensor, resulting in a (k - 1)th order tensor w. Take the best rank-one tensor approximating w. Find a vector t so that t tensor w is the best possible approximation to v, by considering the components of t individually.

Again a reduction, this time to the same problem for k - 1 (which we solve by induction with SVD as the base case) and to finding a scalar t_i  such that t_i * w approximates v_i best. This sounds like a problem that is manageable; I think it's probably convex.

However, this is certainly only an approximation algorithm: for the tensor

 tensor  minus  tensor  

every resulting approximation is the zero matrix. I guess I will implement this and see how well it performs for images...

Erik Postma
Maplesoft. 

I think the best option is selecting such a polygon at random. The probability of three diagonals being convergent is zero in that case. However, you are going to get = 126 points of intersection inside your 9-gon (and twice as many outside it), so they are going to be close together.

Let me give you a code example, and let's also find the minimum distance between two intersections. There's some discussion afterwards.

with(Statistics):
with(simplex):

s := Sample(RandomVariable(Normal(0, 1)));
# This returns a procedure that will create samples for us. (This is underdocumented 
# in Maple 14 but we expect to improve that for the next version.)
hull := []:
while nops(hull) < 9 do
  hull := convexhull([op(hull), convert(s(2), list)]);
end do:

hull_plot := plots:-display(plot(hull, style=point, color=black, symbolsize=20), 
  plot([op(hull), hull[1]], color=green)):
hull_plot;
# This shows the nine points that we've found.

diagonalpairs := map(proc(quad)
  local p1, p2, p3, p4;
    p1, p2, p3, p4 := op(quad);
    return [[p1, p2], [p3, p4]], [[p1, p3], [p2, p4]], [[p1, p4], [p2, p3]];
  end proc, combinat:-choose(hull, 4)):
# We've now found all pairs of disjoint pairs of points defining the hull (intersecting
# inside the convex hull and outside it).

# The following procedure takes two points and returns a procedure mapping an 
# x-value to the corresponding y-value on the line connecting the two points.
# It will error out if it gets two points with the same x-value - but that happens
# with probability 0.
toProc := proc(p1, p2)
  local a, b, x;
    a := (p2[2] - p1[2])/(p2[1] - p1[1]);
    b := p1[2] - a * p1[1];
    return unapply(a*x + b, x);
end proc;
intersections := map(proc(quad)
    local x, x0, p1p, p2p;
    p1p := toProc(op(quad[1]));
    p2p := toProc(op(quad[2]));
    x0 := solve(p1p(x) = p2p(x), x);
    return [x0, p1p(x0)];
  end proc, diagonalpairs):
# We now have the 378 points of intersection of each pair of diagonals.

# Display the hull and the intersections (viewport shrunk to fit the hull).
plots:-display(hull_plot, plot(intersections, style=point), 
  view=(min .. max)~(map2(map2, op, [1, 2], hull)));

mindist := sqrt(min(map(x -> (x[1][1] - x[2][1])^2 + (x[1][2] - x[2][2])^2, 
  combinat:-choose(intersections, 2))));

The minimum distance was approximately 0.0024 in the case that I tried.

An interesting question (to me) is what a good distribution is to choose the points from, if you want to obtain a convex hull consisting of 9 points by selecting as few points as possible (in expectation). It seems likely that somebody has thought about this before. Does anyone know of a paper? (Update: you can get some very interesting distributions of the points over the interior if you select a uniform distribution, because the polygon becomes an approximation to a square.)

Finally: maybe you were hoping to get nice rationals as a result. You could of course just use the floating point numbers as rationals, but that is not exactly nice. Another option is to try rounding the coordinates of the points in the hull using convert(..., confrac); that is, taking one of the convergents in the continued fraction expansion. Then of course you should check afterwards that the minimum distance is still positive: too aggressive rounding will make the intersections be equal. Let us know if you need help with this.

Hope this helps,

Erik Postma
Maplesoft.

Hi Glen,

One useful source of information on data structures is ?MaplePortal/Tutorial6. The link points to the Maple 14 version of that help page; I don't think it existed in Maple 12, but most of the information there also holds for that version.

Another starting point is the ?Array help page, if you do decide on using that data structure - which, I agree with Joe, sounds like a good idea from your description. The third help page I would recommend is ?rtable_indexing. If you're very concerned about performance, then generally your reading list should include some of the ArrayTools help pages, such as the ones for ?ArrayTools/Alias, ?ArrayTools/Copy, and ?ArrayTools/BlockCopy, but for you maybe ?rtable_scanblock would be even more useful.

The other mutable data structure is the ?table (essentially a hash table), but that is more suitable if you need to store information indexed by things not easily expressible as tuples of integers.

As for your remark about using for versus while: in Maple you don't need to choose - you can use both in the same command! See ?do.

Hope this helps,

Erik Postma
Maplesoft.

Hi sund002,

It looks like the last line there is not considered part of the input of the Maple Document you've created. Not having a very good understanding of how the GUI works, and also not knowing how you created the line of input that it is on, I can't tell you exactly what's going on, but you can get around this issue by copying that line of input and pasting it into a new "execution group" using the "[>" button or the shortcut keys Control-J and Control-K (on Linux and, I believe, Windows).

If you like to use the "red Maple" or "Maple Input" style of input, then you might want to create worksheets instead of documents in the future. These have fewer features but are somewhat more predictable. You can try this out by going to File -> New -> Worksheet Mode and make it the default by going to Tools -> Options... -> Interface -> Default format for new worksheets -> select Worksheet -> Apply Globally. Under those same conditions (liking to work with "red Maple"), you might want to select Tools -> Options... -> Display -> Input Display -> select Maple Notation -> Apply Globally, if you haven't done that yet.

Hope this helps,

Erik Postma
Maplesoft.

Hi RNiven,

If I do

plotsetup(eps, plotoutput="temp.eps", plotoptions="axiswidth=10,axisheight=500");      
plot(x,x=0..1);

I get a plot that is very skewed, so certainly something is affected. This is on linux, not on mac, admittedly. ImageMagick identifies the full image, including margins, as being 624x38; I ran identify temp.eps from the linux command line to obtain that info. Specifying other sizes leads to numbers that seem correct as well. What is the exact command you are running that do not seem to have an effect?

I think that dual axis plots are not supported on all plot devices; this is documented on the ?plots/dualaxisplot page. (The last line of the Description section.) I'm afraid you don't have a lot of options there, sorry.

Hope this helps,

Erik Postma
Maplesoft.

Great answers there already. I'd just like to add one more approach that may not be the most useful here, but that I've found to work very well for other cases; especially if you're dealing with potentially big matrices and you (therefore) care about efficiency. This approach is use of the ?ArrayTools/Alias command. You could use it here as follows:

TransposeByChangingOrdering := proc(m :: Matrix, $)
  return ArrayTools:-Alias(m, ListTools:-Reverse([LinearAlgebra:-Dimensions(m)]),
    `if`(rtable_options(m, 'order') = 'Fortran_order', 'C_order', 'Fortran_order'));
end proc;
A := TransposeByChangingOrdering(Matrix(A));

Note that performing linear algebra on C order matrices is typically much slower than on Fortran order matrices, because Maple needs to convert to Fortran order first in order to enable use of the fast BLAS routines. This is only relevant for datatype = float matrices, though, I think. (You can always check what's going on by setting infolevel[LinearAlgebra] to a suitably high value.)

Hope this helps,

Erik Postma
Maplesoft.

Hi widmer,

For your first question, I would suggest the use of the ?assuming facility: you can tell Maple that you want the answer for real (as opposed to complex nonreal) values of t, as follows:

PrincipalNormal(r(t), t, normalized) assuming t :: real;

You can turn off the "basis format" display (as it's called) using ?VectorCalculus/BasisFormat, as follows:

with(VectorCalculus):
BasisFormat(false);

After running the BasisFormat command, your vectors will display in the "normal" format. If you've loaded the VectorCalculus package already (as suggested by your use of PrincipalNormal without the VectorCalculus prefix), then of course you don't have to load it again.

Hope this helps,
Erik Postma.

Hi tropic,

You will first want to store the data together in a list or ?Matrix. For example, like this, if you want to keep the loop you had:

data := Matrix(0, 0);
for i to 4 do
  x := 0.01*i;
  if (i = 1) then
    y := 0.2;
    z := 0.05:
  elif (i = 2) then
    y := 0.3;
    z := 0.5:
  elif (i = 3) then
    y := 0.5;
    z := 1.05:
  elif (i = 4) then
    y := 1.2;
    z := 2.05:
  else
  end if;

  mm(i, 1) := x;
  mm(i, 2) := y;
  mm(i, 3) := z;
end do;

Or even easier for this case:

mm := Matrix([[0.01, 0.02, 0.03, 0.04], [0.2, 0.3, 0.5, 1.2], [0.05, 0.5, 1.05, 2.05]]);

In both cases, you can make the plots as follows:

p1 := plots:-pointplot(mm[1], mm[2], style=line, color=blue, legend=y):
p2 := plots:-pointplot(mm[1], mm[3], style=line, color=red, legend=z):
plots:-display(p1, p2, axes=boxed, scaling=constrained, labels=[x, ""]);

I changed style=point to style=line, because you said you wanted the plot drawn as two lines. And I would actually leave out the scaling=constrained option from this plot, because it makes it very hard to see what's going on. If you do indeed leave that out, you get this plot:

Hope this helps,

Erik Postma
Maplesoft.

We might be able to help you better if you give us the exact equations you are trying to solve.

Other than that, there is only one issue I can think of, given what you wrote: in Maple, xy is a single variable that happens to have a name consisting of two characters; for the multiplication of x and y, you can explicitly enter x * y or, in 2D input mode, enter a space in between, as x y. The exponents should work the way you entered them here.

Erik Postma
Maplesoft.

Have you looked at the ?plot/tickmarks help page? Using the "list of equations" form, you can control the appearance of the tick marks in a pretty precise manner. This does mean that you have to decide where the tick marks are going to be placed.

Hope this helps,

Erik Postma
Maplesoft.

Hi!

The issue is that you instruct Maple to compute the minimum of all 1000 entries in every iteration of the loop, even at the beginning, when most of the values y[k] do not have a value yet. If you do this instead:

for i from 1 to 1000 do
  x := 1 + i;
  y[i] := 2*x + 1;
od;
a := min(seq(y[k], k=1..1000));
print(a);

then you should get the answer (5).

For this case, it is going to be even more efficient, and even simpler, if you write the values of y[k] inline in the seq statement, without the loop beforehand. So you would write

a := min(seq(2 * (1 + i) + 1, i = 1 .. 1000));

instead of the whole loop above. But if the computation of y[k] is more complicated and involves, for example, values y[m] for some smaller values m < k, then the scheme with the loop could be easier.

Hope this helps,

Erik Postma
Maplesoft.

Try this (with thanks to John's baby names document) :

    q := cat("GET /SizedImages/image.ashx?file=991a3b8f2151d4732c78ce6920b9dc88_60_60.jpg HTTP/1.0\n",
             "Host:    www.mapleprimes.com\n\n");
    sock := Sockets:-Open("www.mapleprimes.com", 80);
    
    Sockets:-Write(sock, q);
    
    s := Array(1..100, datatype=integer[1]);
    p := Sockets:-ReadBinary(s, sock):
    n := p;
    while p = 100 do
        a := Array(1..100, datatype=integer[1]);
        p := Sockets:-ReadBinary(a, sock):
        s(n + 1 .. n + p) := a[1 .. p];
        n := n + p;
    end do:
    
    Sockets:-Close(sock);
 
    k := 1;
    while s[k] <> 13 or s[k+1] <> 10 or s[k+2] <> 13 or s[k+3] <> 10 do
      k := k + 1;
    end do:

    s := s[k+4 ..]:

    FileTools:-Binary:-Write("file.jpg", s);
    FileTools:-Binary:-Close("file.jpg");

    img := ImageTools:-Read("file.jpg");
    ImageTools:-View(img);

Hope this helps,

Erik Postma
Maplesoft.

PieChart and the like should really have a legend option, and maybe I will at some point get around to that. Until then, here is a command that will match the sectors up with the left hand sides of your list T:

Q := subs(zip((pol, legend) -> pol = op(0, pol)(op(pol), LEGEND(legend)), select(type, [op(P)], specfunc(anything, POLYGONS)), map(lhs, T)), P):
plots[display](Q, title = ["Distribución de las exportaciones de Chile en el año 2009", font = [TIMES, BOLD, 18]], font = [TIMES, BOLD, 12]);

This is not necessarily guaranteed to work in all cases, I'm afraid. Making a legend for the PieChart plot structure after the plot has been generated cannot really be done reliably. But still, it will work in a large number of cases in versions of Maple that are currently available - and it works for this particular plot in Maple 12.

Hope this helps,

Erik Postma
Maplesoft.