epostma

1364 Reputation

17 Badges

12 years, 305 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 replies submitted by epostma

+1, nice solution. I must admit I didn't try to optimize for speed at all since the demo input [2, 1, 4] was so small - but that's not likely to be the problem that Markiyan wanted to solve :)

It would probably be possible to generate a Compiler:-Compile'able version with some effort, if the highest possible speed is required (iterating over the integers 0 .. 2^n - 1 and using the binary representation to give the signs of the coefficients - with a maximum of 63 or so for n, but that will take way too long anyway), but it's probably better to invest the effort in making a good branch-and-bound strategy. The problem is NP-complete, since Knapsack reduces to it, so solving it perfectly for truly large problem sizes is not going to work.

Erik Postma
Maplesoft.

Indeed the condition m=0 is the problem here: Maple is correctly telling us that there is, in general, no solution to that combination of equations. If we substitute m = 0 in both equations, then we get the solution x(t) = y2(t) = 0 for all t. Otherwise, if we leave m symbolic as is, we can get a solution if we leave the initial condition x(0) = 0 out; it's kind of complicated so I won't include it here. Leaving out the condition y2(0) = 0 still leaves an inconsistent system.

Hope this helps,

Erik Postma
Maplesoft.

Indeed the condition m=0 is the problem here: Maple is correctly telling us that there is, in general, no solution to that combination of equations. If we substitute m = 0 in both equations, then we get the solution x(t) = y2(t) = 0 for all t. Otherwise, if we leave m symbolic as is, we can get a solution if we leave the initial condition x(0) = 0 out; it's kind of complicated so I won't include it here. Leaving out the condition y2(0) = 0 still leaves an inconsistent system.

Hope this helps,

Erik Postma
Maplesoft.

A mini addition to your nice answer: I realized as I saw your code that you can now get some of the speed up of in-kernel code for filling Vectors with a constant, without using ?ArrayTools/Fill, as follows:

(**) bl := Vector(10^7, datatype=float[8]):
(**) kernelopts(bytesused);                 
                                                        82026976

(**) bl[1..5*10^6] := 1:
(**) kernelopts(bytesused);
                                                        82041496

(**) bl[..] := 2:
(**) kernelopts(bytesused);
                                                        82055672

So if you assign a constant to a subrange of an rtable, then no additional allocations are done. (ArrayTools:-Fill is still faster, though, presumably because it is more specialized: you could use the assignment statement for more complicated subsets of the data than ArrayTools:-Fill.)

Hope this helps,

Erik Postma
Maplesoft.

A mini addition to your nice answer: I realized as I saw your code that you can now get some of the speed up of in-kernel code for filling Vectors with a constant, without using ?ArrayTools/Fill, as follows:

(**) bl := Vector(10^7, datatype=float[8]):
(**) kernelopts(bytesused);                 
                                                        82026976

(**) bl[1..5*10^6] := 1:
(**) kernelopts(bytesused);
                                                        82041496

(**) bl[..] := 2:
(**) kernelopts(bytesused);
                                                        82055672

So if you assign a constant to a subrange of an rtable, then no additional allocations are done. (ArrayTools:-Fill is still faster, though, presumably because it is more specialized: you could use the assignment statement for more complicated subsets of the data than ArrayTools:-Fill.)

Hope this helps,

Erik Postma
Maplesoft.

For what it's worth, I've entered this bug into our internal tracking system; that's not a promise that I will have time to do something about it soon, but I'll do my best.

Erik Postma
Maplesoft.

For what it's worth, I've entered this bug into our internal tracking system; that's not a promise that I will have time to do something about it soon, but I'll do my best.

Erik Postma
Maplesoft.

One more little tweak: starting from acer's last version, replace the line

data[y, ix] := data[y, ix] + 1;

by

data[y, ix] := data[y, ix] + iy^(-0.5);

This emphasizes the first couple of iterations over the later ones, and thus shows the 'trails' better as the expression converges towards the attractor. See below. Of course iy^(-0.5) can be replaced by any other function of iy; maybe you'd want to emphasize the first couple of iterations even more, or you'd want to emphasize a few iterations around iy=500, or...

Erik.

@acer Awesome :) I'd like to lighten up the darker areas a bit; so I inserted

map[evalhf,inplace](x -> x^0.8, img_h); 

immediately after the call to worker. Using 50000 iterations per column instead of 500 (using a measly 8 seconds) I got this:

@acer : how is this?

with(ImageTools):
BF8 := proc(ra, rb, xr, yr, n)
local data;
  data := ImageTools:-Create(yr, xr, 1);
  worker(data, HFloat(evalf(ra)), HFloat(evalf(rb)), xr, yr, n);
  return data;
end proc:
worker := proc(
  data :: Array(datatype=float[8], order=C_order),
  ra :: float[8],
  rb :: float[8],
  xr :: posint,
  yr :: posint,
  n :: posint)
local
  ix :: posint,
  iy :: posint,
  r :: float[8],
  rscale :: float[8],
  x :: float[8],
  y :: posint;

  rscale := (rb - ra) / (xr - 1);
  for ix to xr do
    r := ra + rscale * (ix - 1);
    x := HFloat(0.5);
    for iy to 3 do
      # Warmup.
      x := r * x * (1-x);
    end do;
    for iy to n do
      x := r * x * (1-x);
      y := trunc(x * yr) + 1;
      data[y, ix] := data[y, ix] + 1;
    end do;
  end do;
end proc:
worker := Compiler:-Compile(worker):
m := BF8(2.6, 4, 1000, 1000, 250):
map[inplace, evalhf](`/`, m, max(m)):
m3 := m ^~ (1/3):
ImageTools:-Write("abc.tiff", m3):

Leads to this image (converted to png). As you see it's upside down; easy to fix but I didn't bother. If your screen is too narrow (as mine is - I use my monitor in portrait mode), there's a scaled down version below.

Erik.

@Joe Riel : I think I like the idea of a 'let ... in ...' statement and I'd probably use it frequently. It would potentially lead to very high levels of nesting, I imagine, which would be an issue for aesthetical reasons - I would probably request to have your maplev-mode not indent the blocks inside these statements. I wonder whether the memory manager could be made to take advantage of the immutability of these guys. Maybe it would also be nice to have a way to close multiple let... in... scopes at once. (It's interesting to note the difference with the ?use statement, by the way, which creates only a temporary binding (not even a real name) for an expression - it doesn't evaluate the expression at the beginning of the statement, but does so every time it is referenced.)

On the whole I agree with Jacques' comments regarding option remember. I think we now have a slightly better approach with ?option/cache in that it can't grow arbitrarily out of control as easily in terms of memory usage overhead, but it can still be tricky in terms of memory access overhead (and other overhead, and considerations such as mutable data structures and side effects).

Erik.

Joe: I think the reason that unapply(f(x), x) doesn't return f is, that the result of unapply(f(x), x) should require a first argument and ignore the second and further arguments. f itself may not have that property. For example:

f := proc(a, b, c, $)
  if not type([_passed], 'list(numeric)') then
    return 'procname'(_passed);
  elif _npassed = 0 then
    return 1;
  elif _npassed = 1 then
    return a;
  else
    return `*`(_passed);
  end if;
 end proc;

Now f(x) evaluates to itself; and unapply(f(x), x) evaluates to, effectively, the identity function on a single argument. However, f itself is a more complicated beast, and maybe unapply was used specifically to 'tame' it.

Erik Postma
Maplesoft. 

Joe: I think the reason that unapply(f(x), x) doesn't return f is, that the result of unapply(f(x), x) should require a first argument and ignore the second and further arguments. f itself may not have that property. For example:

f := proc(a, b, c, $)
  if not type([_passed], 'list(numeric)') then
    return 'procname'(_passed);
  elif _npassed = 0 then
    return 1;
  elif _npassed = 1 then
    return a;
  else
    return `*`(_passed);
  end if;
 end proc;

Now f(x) evaluates to itself; and unapply(f(x), x) evaluates to, effectively, the identity function on a single argument. However, f itself is a more complicated beast, and maybe unapply was used specifically to 'tame' it.

Erik Postma
Maplesoft. 

What is the exact message you get? Is it "Error, recursive assignment"? If so, then you did not execute the statement c := 0, or you assigned something else to c after assigning 0 to c. Executing those two statements in order cannot give you that error message.

Hope this helps,

Erik Postma
Maplesoft.

@acer : Yes, I meant to include that option. I wonder where I dropped it - I ran the code, so I must have done something wrong copying and pasting. I agree with your other remarks as well.

Erik.

2 3 4 5 6 7 8 Last Page 4 of 20