epostma

1394 Reputation

17 Badges

13 years, 134 days
Maplesoft

Social Networks and Content at Maplesoft.com

Maple Application Center
I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity


These are answers submitted by epostma

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 <9,4> = 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:

PLOT(CURVES(Matrix(4, 2, {(1, 1) = HFloat(0.100000000000000002e-1), (1, 2) = HFloat(.200000000000000011), (2, 1) = HFloat(0.200000000000000004e-1), (2, 2) = HFloat(.299999999999999989), (3, 1) = HFloat(0.299999999999999989e-1), (3, 2) = HFloat(.500000000000000000), (4, 1) = HFloat(0.400000000000000008e-1), (4, 2) = HFloat(1.19999999999999996)}, datatype = float[8], storage = rectangular, order = Fortran_order, shape = []), LEGEND(y), COLOUR(RGB, 0., 0., 1.00000000), STYLE(LINE)), CURVES(Matrix(4, 2, {(1, 1) = HFloat(0.100000000000000002e-1), (1, 2) = HFloat(0.500000000000000028e-1), (2, 1) = HFloat(0.200000000000000004e-1), (2, 2) = HFloat(.500000000000000000), (3, 1) = HFloat(0.299999999999999989e-1), (3, 2) = HFloat(1.05000000000000004), (4, 1) = HFloat(0.400000000000000008e-1), (4, 2) = HFloat(2.04999999999999982)}, datatype = float[8], storage = rectangular, order = Fortran_order, shape = []), LEGEND(z), COLOUR(RGB, 1.00000000, 0., 0.), STYLE(LINE)), AXESLABELS(x, ""), AXESSTYLE(BOX))

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.

Hi herclau,

I'm not sure what factor exactly you want to factor out - this one example does not provide enough information to define that uniquely. What would intuitively make sense to me would be to factor out the GCD of the entries of the vector, as follows:

myFactorization1 := proc(v :: rtable, $)
local divisor;
  divisor := foldr(gcd, op(convert(v, list)));
  return divisor, v / divisor;
end proc;

but for the case above, that would return b, <0, 1, 0, 5>. Maybe you would want the LCM of the nonzero entries of the vector? That could be done as follows:

myFactorization2 := proc(v :: rtable, $)
local divisor;
  divisor := foldr(lcm, op(remove(`=`, convert(v, list), 0)));
  return divisor, v / divisor;
end proc;

If you expect to find vectors from a zero-dimensional vector space then you need to adapt both routines, and if you expect all-zero vectors then you need to adapt the last one of these two for that.

There's no really nice solution for getting this as a true factorization that is a product of things, because multiplying a vector by a scalar always evaluates the product. You could instead have the procedures return `%*`(divisor, v / divisor). This is an expression that will turn into the original product once you apply 'value' to it:

a := `%*`(b, c); # this will remain unevaluated through all computations that don't apply value
value(1 + a*d); # this returns 1 + b*c*d

Hope this helps,

Erik Postma
Maplesoft.

Hi Neill,

You've already found that you can do rhs(Eq1) = lhs(Eq1). Just a tiny bit shorter is this:

algsubs((rhs = lhs)(Eq1), Eq2);

This works because if you use something like f = g as a function applied to x, then this evaluates to f(x) = g(x).

As you already suggest in the title, you can also use algsubs. Whenever the thing you are replacing is not just a variable name, algsubs is more reliable than subs. So I'd recommend using that for this case. An alternative would be to roll your own version of algsubs using solve or isolate.

Hope this helps,

Erik Postma
Maplesoft.

Hi afeddersen,

Sorry - I seem to have a habit of answering direct questions with "here is something different you could do that might accomplish your goals in a better manner", and I'm afraid this is going to be such an answer.

The label feature in Maple is great, but it mostly works with direct user interaction: there's no easy access to programmatically retrieve the value of label (x.y) for example, if x and y are integers; instead you have to enter these numbers manually. Loops, of course, are very much programmatic constructs, and thus these two don't mesh very well.

One option would be to just store whatever value is printed out in a ?table or an ?rtable (such as a ?Vector). So if you have a loop that looks like this:

running_total := 0;
for i to 25 do
  running_total := running_total + i:
  i^2 - running_total;
end do;

to print the numbers i*(i-1)/2, you could instead decide to store these in a Vector:

results := Vector(25);
running_total := 0;
for i to 25 do
  running_total := running_total + i:
  results[i] := i^2 - running_total;
end do;

Now you can access the ith value as results[i] - with the added advantage over labels that it also works programmatically, not just when you type it directly. Does that help for what you're trying to do?

Erik Postma
Maplesoft.

4 5 6 7 8 9 10 Page 6 of 11