@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),
  %
               );

Fantastic, thanks.

acer

Fantastic, thanks.

acer

@thwle Robert's code uses a syntax designed for the command plots:-arrow

There is another command, plottools:-arrow, and when you load `plottools` after loading `plots` (using `with`) then you are clobbering the previous binding of `arrow`.

You can force it to use the intended `arrow`, by explicitly using its so-called long-form,

with(plots): with(plottools):

animate(plots:-arrow,[[0,0],[0,cos(t)]],t=0..4*Pi,view=[-1..1,-1..1]);

@thwle Robert's code uses a syntax designed for the command plots:-arrow

There is another command, plottools:-arrow, and when you load `plottools` after loading `plots` (using `with`) then you are clobbering the previous binding of `arrow`.

You can force it to use the intended `arrow`, by explicitly using its so-called long-form,

with(plots): with(plottools):

animate(plots:-arrow,[[0,0],[0,cos(t)]],t=0..4*Pi,view=[-1..1,-1..1]);

@Robert Israel I noticed that a single plot, as the kludge insertion, doesn't include the text portion which shows the current value of the animation paraneter. Ie, a textplot with "z=...." in the plot region.

If you can discern that absence, for the kludged insertion, then you might be able to instead make the kludge be a 1- (or unchanging 2-frame) animation that also computes and inserts quickly.