@carizard Try spelling it correctly, as VectorCalculus instead of VektorCalculus.

@Bruce Scranton 

It sounds as if your particular use case is not so much to get a very specific pixel-size for your plot as it is to just get it inserted by the GUI as large. And you have precomputed the plot structure, and even have it saved to a name.

If that's so, then this might help as an alternative approach. The Standard GUI can remember the size of a plot output region, and re-use that size when overwriting with a recalculated plot. I believe that this may even be turned on/off with the option checkbox labeled "Remember plot attributes when re-executing", in Tools->Options->Interface from the GUI's main menubar. You can enable this, and try it instead to get what you're after.

So insert a new execution group (in your Worksheet) or a new paragraph (in your Document). Plot something simple, like `plot(sin(x),x=...)`, in that execution group. Resize the resulting displayed plot output region, using the mouse-pointer. Then copy-and-paste something else into the input line, overwriting the `plot(sin(x),x=...)` call. Replace it by pasting in something else that would normally result in a 2D plot being displayed, such as kochstar(i) or plots:-display(kochstar(i)). Then execute that input line. The result should "re-use" the same plot output region, which you have already re-sized as desired.

acer

@PatrickT Here is one (crude) way, using subs. There are no offsets to make handling the indices nicer, though, so it doesn't give much more functionality.

restart:

BF := proc(fexpr,ra,rb,xr,yr,n,init::float)
local img_h,img_s,img_v,worker;
   uses LinearAlgebra, ImageTools, Compiler;
   worker := subs(':-DUMMY'=fexpr,
                  proc(data::Array(datatype=float[8],order=Fortran_order),
                       ra::float[8],rb::float[8],xr::posint,yr::posint,
                       n::posint,init::float[8])
                  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 := init;
                     for iy to 0 do
                        x := DUMMY;
                     end do;
                     for iy to n do
                        x := DUMMY;
                        y := trunc(x * yr) + 1;
                        if ix>=0 and xr>=ix and y>=0 and yr>=y then
                           data[y, ix] := data[y, ix] + 1;
                        end if;
                     end do;
                  end do;
                  NULL;
               end proc);
   worker := Compile(worker):
   img_h := Array(1..xr,1..yr,'datatype'='float[8]','order'='Fortran_order');
   try
      worker(img_h,HFloat(evalf(ra)),HFloat(evalf(rb)),xr,yr,n,init);
   catch:
   end try;
   rtable_options(img_h,'subtype'=Matrix);
   MatrixScalarMultiply(img_h,360.0/max(img_h),'inplace');
   img_v:=Array(img_h);
   rtable_options(img_v,'subtype'=Matrix);
   MatrixScalarMultiply(img_v,1/3.6,'inplace');
   rtable_options(img_h,'subtype'=Array);
   rtable_options(img_v,'subtype'=Array);
   img_s:=Array(1..xr,1..yr,fill=1,'datatype'='float[8]',
                'order'='Fortran_order'):
   HSVtoRGB(CombineLayers(img_h, img_s, img_v)):
end proc:

m := CodeTools:-Usage( BF(r*x*(1-x), 0, 4, 900, 900, 5000, 0.5) ):

ImageTools:-View(m);

@PatrickT 

If the procedure `f` which implements the formula is itself evalhf'able then you can pass it to evalhf(`worker`(..,`f`,..)) and have it be interpreted as fast as evalhf can do, without penalty of escaping evalhf for usual, slower interpreting. But we want `worker` Compiled instead, so that it executes 15 or more times faster even than evalhf(worker(...)).

It can likely be done by either using `subs` to overwrite a version of `worker` that has a DUMMY name as placeholder. Or possibly by taking the expression formula and unapplying it as an `option inline` operator which then becomes fully explicit as a replacement of a function call inside `worker`. (Sorry, that last sentence is hideous.) AFAIK the Compiler does not currently support references from within a given proc out to other Compiled procs (or at least, not with kludging).

More later -- have to rush off now.

@Bruce Scranton 

As Joe mentions, in general you could create a catch-all folder for such add-on material. The key thing there is that certain locations for folders named "lib" will automatically get added to `libname` upon Maple startup.

But the underlying mechanism (read: trick) at work for this sized plot insertion is that it has to first save to an .hdb help file. Offhand, I can't recall whether the private .hdb file it's using for that is referenced with a hard-coded location. I suspect that it tries to write to a particularly named .hdb file, the specific location of which may or may not matter. It's probably safe to try moving it to cat(kernelopts(mapledir),"/toolbox"), as Joe described.

I keep meaning to come back to this thing. I now believe that I might be able to wrap the two relevent INTERFACE_XYZBLAH() calls inside streamcall() and by doing so get the full effect with just a single procedural top-level return. By which I mean I should be able to get rid of that second function call, to insert into the current worksheet.

Ideally the plot should be placed into an output area, when inserted into the running worksheet. And it might be nice to (optionally) allow the actual plot command in question to be inserted into the corresponding input area. Hopefully, that would mean (behaving as usual) that subsequent re-execution of that inserted plot command would re-fill the presized output and retain the specified size. And the `Resize` call could be deleted by the worksheet author, immediately following insertion. The `Resize` call isn't wanted, longterm, as re-execution of the worksheet would cause repeated insertions, not overwrites.

Lastly, food for thought in this holiday season: the underlying mechanism allows on-the-fly insertion of programatically constructed portions of worksheet XML. What other kinds of compound and complex worksheet objects do we want to be able to programatically construct and insert, live, into a running Worksheet or Document? (I have a running wishlist of my own,...)

@Joe Riel 

It works with both the right-click context-menu .eps export and the programmatic `plotsetup` .eps driver when using the Standard GUI.

I did mention the Standard GUI in the first sentence of the Post, but I apologize to all for not explicitly stating that this was not intended to be related to the Classic GUI.

@headbam 

Your conclusion is correct.

See Atomic Identifiers for an alternative was to subscript names so that they are not table member references (and so are fully distinct from the base name).

@headbam 

Your conclusion is correct.

See Atomic Identifiers for an alternative was to subscript names so that they are not table member references (and so are fully distinct from the base name).

@PatrickT 

1) I suspect that the warmup is a attempt to not get a line of zeros as false partial results.

2) A common interpretation is that, on the "color wheel", hue is cyclic. See google for illustration. The ImageTools package interprets hue in this way, at least, in terms of degrees round a circle. Now, the logistic map was incrementing entries of Array `img_h` with each placement. So `img_h` comes out of the `worker` with nonnegint values. Scaling them all down by the maximal value gets them all to the range 0 to 1. Then scaling up by 360 places them around the circle, as hue values for ImageTools. I merged those together, scaling by 360/max(img_h). The scaling of one of the other hsv layer Arrays was just rough, for convenience and taste.

3) LinearAlgebra:-MatrixScalarMultiply calls an external compiled BLAS function, which will be faster than evalhf. But LinearAlgebra currently only accepts Matrices and Vectors and not Arrays. And ImageTools wants Arrays. Fortunately, the `subtype` of an rtable can be toggled, inplace and with no data copying cost, if the indexing matches. You should be able to use this, instead, if you really wanted to not call LinearAlgebra there,

map[evalhf,inplace]( `*`, img_h, 360.0/evalf(max(img_h)) );

I really love how Erik's color image above has sharp, consistent changes in color whenever it bifurcates. But there are lots of choices, such as non-constant scalings for example. It might be nice to analyse the distribution of the raw img_h values, so as to smear them more evenly (once) around the hue wheel.

The code is pretty much hard-coded for a particular map. So one nice improvement would be to rework it to allow the formula as a parameter. That would necessitate some changes to keep the computed discrete values within bounds of the Arrays.

There's a lot of scope for playing around with the colouring.

restart:

BF8 := proc(ra,rb,xr,yr,n)
local img_h,img_s,img_v,img_hsv;
  img_h := Array(1..xr,1..yr,'datatype'='float[8]','order'='Fortran_order');
  worker(img_h,HFloat(evalf(ra)), HFloat(evalf(rb)),xr,yr,n);
  rtable_options(img_h,'subtype'=Matrix);
  LinearAlgebra:-MatrixScalarMultiply(img_h,360.0/max(img_h),'inplace');
  img_v:=Array(img_h);
  rtable_options(img_v,'subtype'=Matrix);
  LinearAlgebra:-MatrixScalarMultiply(img_v,1/3.6,'inplace');
  rtable_options(img_h,'subtype'=Array);
  rtable_options(img_v,'subtype'=Array);
  img_s:=Array(1..xr,1..yr,fill=1,'datatype'='float[8]',
               'order'='Fortran_order'):
  img_hsv:=ImageTools:-CombineLayers(img_h, img_s, img_v):
  return ImageTools:-HSVtoRGB(img_hsv):
end proc:

worker := proc( data::Array(datatype=float[8],order=Fortran_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;
  NULL;
end proc:

worker := Compiler:-Compile(worker):

m := CodeTools:-Usage( BF8(2.6, 4, 1200, 1200, 500) ):

#ImageTools:-Preview(m);
#ImageTools:-View(m);

ImageTools:-Write(cat(kernelopts(homedir),"/bifucolor.jpg"),
                  m,'format'='JPEG'):

 

 

acer

@epostma It is beautiful.

Your code produces the grayscale tiff in under a second on my machine, with very modest memory use. This should be taken as target performance for the Standard GUI, for `plot`+`style=point` with an mx2 float[8] Array as well!

Now, it looks as if you are keeping track of the number of times that you "end up" at any data[y,ix] point, and using that for the grayscale. How about using those values (scaled) for the Hue layer of an HSV image?

If only we could programatically embed images right into worksheets, instead of fumbling with resource hungry point-plots. And of course I don't count ImageTools:-Preview which is a rough scaled-down, slow point-plot, or ImageTools:-View which is only a popup.

acer

I have branched this off as a separate Post. It was originally an Answer to a Question whose author had also asked a duplicate. The source is a modification of from this earlier post. Another related Question did not yet have a satisfying answer (in the efficiency sense, which was mentioned by the submitter).

Patrick T. has started us off. I think that computing directly into a float[8] Array, and then writing that out directly using ImageTools, should be a way to get sharply detailed results quickly. The crucial parts might be Compiled, or possibly Threaded with evalhf, and act inplace on a given Array. But until someone codes that up, here is a modification of Patrick's tweak.

This runs in Maple 15. It also ran in Maple 12, but without the Threads:-Sleep call, whose purpose is only to incur a delay while the plot file driver writes out the first, large image file. In Maple 12, one would have to manually run it, and wait a little. On a fast i7 in Windows 7, 10 seconds seems long enough for the 4096x4096 jpeg.

Bifurcation := proc(initialpoint,xexpr,ra,rb,acc)
  local p1,hr,A,L1,i,j,phi:
  global r,L2:
  hr := unapply(xexpr,x);
  A := Vector(600):
  L1 := Vector(acc*500):
  for j from 1 to acc+1 do
    r := (ra + (j-1)*(rb-ra)/acc):
    A[1] := hr(initialpoint):
    for i from 2 to 500 do
      A[i] := evalf(hr(A[i-1])):
    end do:
    for i from 1 to 400 do
      L1[i+400*(j-1)] := [r,A[i+100]]:
    end do:
  end do:
  L2 := {seq(L1[i], i = 1..acc*400)}:
  p1 := plots:-pointplot(L2, 'symbol' = solidcircle, 'symbolsize' = 2, 'color' = blue):
  unassign('r'):
  return(p1):
end proc:

P1 := Bifurcation(1/2,r*x*(1-x),2.5,4,250):

P:=plots:-display(P1, 'axes' = box, 'labels' = [r, x] ):

# A very large image is needed, to get symbolsize=2 to be seen.
plotsetup(jpeg,plotoptions="height=4096,width=4096",
          plotoutput=cat(kernelopts(homedir),"/bifu4096.jpg")):

plots:-display(P);

Threads:-Sleep(10): # or be patient, and wait

image:=ImageTools:-Read(cat(kernelopts(homedir),"/bifu4096.jpg"),
                        format=JPEG):

ImageTools:-Write(cat(kernelopts(homedir),"/bifu640.jpg"),
                  ImageTools:-Scale(image,1/6.4),format=JPEG):

I'm not altogether content with the above code. I didn't change the methodology. (It builds L2 as a listlist for `pointplot` instead of a float[8] Array for `plot` with style=point. It doesn't use evalhf or any other acceleration. I'm a little picky about globals, in the sense that I prefer code without them, when possible.) It's a good start, to build from, and the topic is interesting.

A few other comments are in order here. I've noticed that very large exported image files are needed for small symbolsize, in order to export something which is non-empty. In fact, the same thing happens in the Standard GUI itself -- if one makes the symbolsize too small, then the plot is shown as non-empty only if the plot window has been (manually) resized as quite large. Is this really useful, helpful, and best?

acer

ps. My firefox is only showing me the Mapleprimes post editor in raw html form, for a day or so. So I don't have the nice editor, and cannot upload files. (Am I the only one with this issue?) If one person who doesn't have this issue would care to produce and upload the final 640x640 image using the above code, that might be nice, thanks.

You don't seem to have shown us exactly how you do want it formatted.

Also, you seem to want x/3 to come out like x/3. or x/.3e1 while at the same time not wanting z/10 to come out as z/.1e2 which I find hard to fathom. Would you accept .33*x and .1*z, or want more precision and if so how much? Or, if you want x/3 and z/10 to be treated differently then how so? How many nonzero trailing decimals can there be in z/K converted to KFLOATINV*z before you'd prefer it as z/KFLOAT?

acer

@herclau 

What Robert has done is construct a transformation procedure, using plottools:-transform, and then apply it to a plot structure. It is procedure application, not multiplication.

Done in steps, it might look like this,

f := x^2-1:
g := -x-1:

origplot:=plot(f-g, x=-1.5 .. 1.5,
               filled=true, color=COLOUR(RGB,.8,.8,.9)):
origplot;

H := unapply([x,y+g],x,y);

transformer := plottools:-transform(H):

transformer(origplot);

plots:-display(
  plot([f,g], x=-1.5 .. 1.5, color=black),
  %
               );

@herclau 

What Robert has done is construct a transformation procedure, using plottools:-transform, and then apply it to a plot structure. It is procedure application, not multiplication.

Done in steps, it might look like this,

f := x^2-1:
g := -x-1:

origplot:=plot(f-g, x=-1.5 .. 1.5,
               filled=true, color=COLOUR(RGB,.8,.8,.9)):
origplot;

H := unapply([x,y+g],x,y);

transformer := plottools:-transform(H):

transformer(origplot);

plots:-display(
  plot([f,g], x=-1.5 .. 1.5, color=black),
  %
               );