acer

32353 Reputation

29 Badges

19 years, 331 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Markiyan Hirnyk The difference lies in inserting Erik's posted 1D Maple Notation code as 2D Math input. it works fine as 1D Maple Notation (plaintext) input. Erik posted 1D Maple Notation code, and you appear to have used it (and reposted it) as 2D Math input.

The following fragment of the posted plaintext 1D code does not work as 2D Math input (if cut and pasted, or if typed manually).

  solvedf := x -> eval(f(:-x), solver(x));

In this particular case, it may be a bug in the 2D Math parser, as even the context-menu cannot convert it from 1D input to 2D input, in situ. This minor variation seems to work, though,

  solvedf := xx -> eval(f(:-x), solver(xx));

More examples of 1D vs 2D parser differences can be seen here.

acer

@Joe Riel Yes, a great deal of it is OS dependent. Windows seems particularly hard hit, here, by the additional entries in libname, in the non-networked case.

But the NFS (or share) mounting is another (larger) additional hit (which may often be less noticed, as it doesn't get reported by the usual `time`of GUI status bar).

What surprised me most about this was that a program which didn't appear to involve i/o in an obvious way could be so affected.

With regard to the OS dependence: years ago, we discussed a Maple Benchmark suite on this site (Mathematica has had an unofficial one, for many years). It never came to pass, but it could be made with a collection of examples which get graded to an aggregate score, so that people can compare the application on various platforms. I'm not sure how subtasks' scores might be weighted.

@LijiH I meant Douglas Harder's Quaternions package. I haven't looked at its source, but I would be surprised if it did not use at least some modern (module, etc) techniques.

I misunderstood what you wanted. Before my first coffee of the day I had guessed that you wanted the special property for any such pair of names you subsequently protected (or somehow tagged as special). I see now that you are likely talking only of the names `a` and `b`.

In that case, you could just make `a` and `b` exports of your module. And then the module's `*` could be coded to treat them specially.

Here is an example that I just ciooked up. It doesn't do anything special with a*a or b*b. It uses a trick for printing: internally the structures it returns are calls to noncommutative `&*`, but since you might not like that operator in printed output it prints using regular infix `*` with names replaced by escaped locals. This printing trick doesn't show up when I cut and paste the output to this post, as plaintext.

restart:

m:=module() option package;
   export `*`, a, b;
   `*`:=proc() #option trace;
   local coeff1, coeff2, coeffprod, x, y;
   if nargs=2 then
      if type(args[1],constant) or type(args[2],constant) then
         return :-`*`(args);
      end if;
      if type(args[1],linear(a)) and type(args[1]/a,constant) then
         coeff1,x:=coeff(args[1],a),a;
      elif type(args[1],linear(b)) and type(args[1]/b,constant) then
         coeff1,x:=coeff(args[1],b),b;
      else
         coeff1,x:=1,args[1];
      end if;
      if type(args[2],linear(a)) and type(args[2]/a,constant) then
         coeff2,y:=coeff(args[2],a),a;
      elif type(args[2],linear(b)) and type(args[2]/b,constant) then
         coeff2,y:=coeff(args[2],b),b;
      else
         coeff2,y:=1,args[2];
      end if;
      coeffprod:=:-`*`(coeff1,coeff2);
      if x=a then
         if y=b then
           return coeffprod;
         elif y=-b then
           return -coeffprod;
         end if;
      elif x=b then
         if y=a then
            return -coeffprod;
         elif x=-a then
            return coeffprod;
         end if;
      end if;
   end if;
   if coeffprod=1 and coeff1=1 then
         return ':-`&*`'(args);
   else
         return :-`&*`(coeffprod,procname(x,y));
   end if;
   return :-`*`(coeffprod,':-`&*`'(op(map(eval,[args]))));
end proc:
end module:

# Could put this in m's ModuleLoad
`print/&*`:=proc()
   :-`*`(op(map(t->`if`(type(t,name),`tools/gensym`(t),t),[args])));
end proc:

# Could put this in m's ModuleLoad
:-`&*`:=m:-`*`:

with(m);
                           [*, a, b]

a*b;
                               1

b*a;
                               -1

a*c*b; # consequence of noncommutativity
                         (a &* c) &* b

a*(c+b);
                          a &* (c + b)

expand(%);
                           1 + a &* c

b*a*b;
                               -b

a*b*a;
                               a

a*(b*a*b);
                               -1

(a*b*a)*b;
                               1

a*b*a*b; # consequence of associativity
                               1

a*4*b;
                               4

a*(4*b);
                               4

b*3*a*6;
                              -18

4*a*b;
                               4

(c+b)*(d+a);
                       (c + b) &* (d + a)

expand(%);
                  c &* d + c &* a - 1 + b &* d

c*a - a*c;
                        c &* a - a &* c

Here's what it looks like in the actual session,

 

restart:

m:=module() option package;

   export `*`, a, b;

   `*`:=proc() #option trace;
   local coeff1, coeff2, coeffprod, x, y;
   if nargs=2 then
      if type(args[1],constant) or type(args[2],constant) then
         return :-`*`(args);
      end if;
      if type(args[1],linear(a)) and type(args[1]/a,constant) then
         coeff1,x:=coeff(args[1],a),a;
      elif type(args[1],linear(b)) and type(args[1]/b,constant) then
         coeff1,x:=coeff(args[1],b),b;
      else
         coeff1,x:=1,args[1];
      end if;
      if type(args[2],linear(a)) and type(args[2]/a,constant) then
         coeff2,y:=coeff(args[2],a),a;
      elif type(args[2],linear(b)) and type(args[2]/b,constant) then
         coeff2,y:=coeff(args[2],b),b;
      else
         coeff2,y:=1,args[2];
      end if;
      coeffprod:=:-`*`(coeff1,coeff2);
      if x=a then
         if y=b then
           return coeffprod;
         elif y=-b then
           return -coeffprod;
         end if;
      elif x=b then
         if y=a then
            return -coeffprod;
         elif x=-a then
            return coeffprod;
         end if;
      end if;
   end if;
   if coeffprod=1 and coeff1=1 then
         return ':-`&*`'(args);
   else
         return :-`&*`(coeffprod,procname(x,y));
   end if;
   return :-`*`(coeffprod,':-`&*`'(op(map(eval,[args]))));
end proc:
end module:

# Could put this in m's ModuleLoad
`print/&*`:=proc()
   :-`*`(op(map(t->`if`(type(t,name),`tools/gensym`(t),t),[args])));
end proc:

# Could put this in m's ModuleLoad
:-`&*`:=m:-`*`:

with(m);

a*b;

b*a;

a*c*b; # consequence of noncommutativity

a*(c+b);

expand(%);

b*a*b;

a*b*a;

a*(b*a*b);

(a*b*a)*b;

a*b*a*b; # consequence of associativity

a*4*b;

a*(4*b);

b*3*a*6;

4*a*b;

(c+b)*(d+a);

expand(%);

c*a - a*c;

 

 

Download starfun.mw

@LijiH I meant Douglas Harder's Quaternions package. I haven't looked at its source, but I would be surprised if it did not use at least some modern (module, etc) techniques.

I misunderstood what you wanted. Before my first coffee of the day I had guessed that you wanted the special property for any such pair of names you subsequently protected (or somehow tagged as special). I see now that you are likely talking only of the names `a` and `b`.

In that case, you could just make `a` and `b` exports of your module. And then the module's `*` could be coded to treat them specially.

Here is an example that I just ciooked up. It doesn't do anything special with a*a or b*b. It uses a trick for printing: internally the structures it returns are calls to noncommutative `&*`, but since you might not like that operator in printed output it prints using regular infix `*` with names replaced by escaped locals. This printing trick doesn't show up when I cut and paste the output to this post, as plaintext.

restart:

m:=module() option package;
   export `*`, a, b;
   `*`:=proc() #option trace;
   local coeff1, coeff2, coeffprod, x, y;
   if nargs=2 then
      if type(args[1],constant) or type(args[2],constant) then
         return :-`*`(args);
      end if;
      if type(args[1],linear(a)) and type(args[1]/a,constant) then
         coeff1,x:=coeff(args[1],a),a;
      elif type(args[1],linear(b)) and type(args[1]/b,constant) then
         coeff1,x:=coeff(args[1],b),b;
      else
         coeff1,x:=1,args[1];
      end if;
      if type(args[2],linear(a)) and type(args[2]/a,constant) then
         coeff2,y:=coeff(args[2],a),a;
      elif type(args[2],linear(b)) and type(args[2]/b,constant) then
         coeff2,y:=coeff(args[2],b),b;
      else
         coeff2,y:=1,args[2];
      end if;
      coeffprod:=:-`*`(coeff1,coeff2);
      if x=a then
         if y=b then
           return coeffprod;
         elif y=-b then
           return -coeffprod;
         end if;
      elif x=b then
         if y=a then
            return -coeffprod;
         elif x=-a then
            return coeffprod;
         end if;
      end if;
   end if;
   if coeffprod=1 and coeff1=1 then
         return ':-`&*`'(args);
   else
         return :-`&*`(coeffprod,procname(x,y));
   end if;
   return :-`*`(coeffprod,':-`&*`'(op(map(eval,[args]))));
end proc:
end module:

# Could put this in m's ModuleLoad
`print/&*`:=proc()
   :-`*`(op(map(t->`if`(type(t,name),`tools/gensym`(t),t),[args])));
end proc:

# Could put this in m's ModuleLoad
:-`&*`:=m:-`*`:

with(m);
                           [*, a, b]

a*b;
                               1

b*a;
                               -1

a*c*b; # consequence of noncommutativity
                         (a &* c) &* b

a*(c+b);
                          a &* (c + b)

expand(%);
                           1 + a &* c

b*a*b;
                               -b

a*b*a;
                               a

a*(b*a*b);
                               -1

(a*b*a)*b;
                               1

a*b*a*b; # consequence of associativity
                               1

a*4*b;
                               4

a*(4*b);
                               4

b*3*a*6;
                              -18

4*a*b;
                               4

(c+b)*(d+a);
                       (c + b) &* (d + a)

expand(%);
                  c &* d + c &* a - 1 + b &* d

c*a - a*c;
                        c &* a - a &* c

Here's what it looks like in the actual session,

 

restart:

m:=module() option package;

   export `*`, a, b;

   `*`:=proc() #option trace;
   local coeff1, coeff2, coeffprod, x, y;
   if nargs=2 then
      if type(args[1],constant) or type(args[2],constant) then
         return :-`*`(args);
      end if;
      if type(args[1],linear(a)) and type(args[1]/a,constant) then
         coeff1,x:=coeff(args[1],a),a;
      elif type(args[1],linear(b)) and type(args[1]/b,constant) then
         coeff1,x:=coeff(args[1],b),b;
      else
         coeff1,x:=1,args[1];
      end if;
      if type(args[2],linear(a)) and type(args[2]/a,constant) then
         coeff2,y:=coeff(args[2],a),a;
      elif type(args[2],linear(b)) and type(args[2]/b,constant) then
         coeff2,y:=coeff(args[2],b),b;
      else
         coeff2,y:=1,args[2];
      end if;
      coeffprod:=:-`*`(coeff1,coeff2);
      if x=a then
         if y=b then
           return coeffprod;
         elif y=-b then
           return -coeffprod;
         end if;
      elif x=b then
         if y=a then
            return -coeffprod;
         elif x=-a then
            return coeffprod;
         end if;
      end if;
   end if;
   if coeffprod=1 and coeff1=1 then
         return ':-`&*`'(args);
   else
         return :-`&*`(coeffprod,procname(x,y));
   end if;
   return :-`*`(coeffprod,':-`&*`'(op(map(eval,[args]))));
end proc:
end module:

# Could put this in m's ModuleLoad
`print/&*`:=proc()
   :-`*`(op(map(t->`if`(type(t,name),`tools/gensym`(t),t),[args])));
end proc:

# Could put this in m's ModuleLoad
:-`&*`:=m:-`*`:

with(m);

a*b;

b*a;

a*c*b; # consequence of noncommutativity

a*(c+b);

expand(%);

b*a*b;

a*b*a;

a*(b*a*b);

(a*b*a)*b;

a*b*a*b; # consequence of associativity

a*4*b;

a*(4*b);

b*3*a*6;

4*a*b;

(c+b)*(d+a);

expand(%);

c*a - a*c;

 

 

Download starfun.mw

It is possible (but a little tricky) to write a Library archive ripper, which runs through procs and modules (and corresponding module members) listed in a .mla file and writes out the source to individual source files.

A large number of routines and packages can be extracted in this way. There are exceptions due to anonymous constructs, etc, but quite a lot of the routines can be extracted to a plaintext form which even be subsequently reloaded successfully.

One can learn a lot about Maple, by reading and studying such source code.

acer

@LijiH 

"1, it will restart the workspace."

No.


"2, it will use packages PolynomialTools and RandomTools."

If you mean that proc1 and proc2 will be codes with calls to those packages, then yes. Just incorporate `uses` statements inside the body of proc1 and proc2. No, you should not try to accomplish this rebinding in a global sense, merely by invoking your package. See the description of the help-page ?proc for a note on `uses`.


"3, it will protect the letters I, J and K."

Yes, sure. Just put calls to the `protect` command inside the ModuleLoad procedure of your package module. You will only see the effects of ModuleLoad being run if you `savelib` the package to a Library archive and then restart and load it.


"4, define M, which is not a proc but will be used in some proc that Helloworld offers."

Sure, a module can have both locals as well as exports. You can define it in the definitionof the module. You can even have proc1 or proc2 write values to such module locals.


You are going to have to create a Library archive, and set libname, to get this to behave how you want. See here for some notes on savelib'ing a module.

@LijiH 

"1, it will restart the workspace."

No.


"2, it will use packages PolynomialTools and RandomTools."

If you mean that proc1 and proc2 will be codes with calls to those packages, then yes. Just incorporate `uses` statements inside the body of proc1 and proc2. No, you should not try to accomplish this rebinding in a global sense, merely by invoking your package. See the description of the help-page ?proc for a note on `uses`.


"3, it will protect the letters I, J and K."

Yes, sure. Just put calls to the `protect` command inside the ModuleLoad procedure of your package module. You will only see the effects of ModuleLoad being run if you `savelib` the package to a Library archive and then restart and load it.


"4, define M, which is not a proc but will be used in some proc that Helloworld offers."

Sure, a module can have both locals as well as exports. You can define it in the definitionof the module. You can even have proc1 or proc2 write values to such module locals.


You are going to have to create a Library archive, and set libname, to get this to behave how you want. See here for some notes on savelib'ing a module.

@TomM You wrote, "However, the whole package is not bundled with Maple, in the sense that the code (dlls) comes with the Maple program distribution." Are you claiming that the dll is shipped with Maple, or not? A dynamic shared CLAPACK library object is, in fact, included in the Maple installation on at least 32 and 64 bit MS-Windows and Linux.

The real and complex double precision funtions (Dxxxxx and Zxxxxx, but not single precision Cxxxxx and Sxxxxx) of CLAPACK 3.0 should be in the clapack.dll that is indeed shipped with Maple. If there are some functions in the dll which are not exported, then that would be something to fix.

There isn't a dynamic BLAS library object shipped with Maple to which you can dynamically link, expect (ATLAS) on Linux.

I agree with you that there is some extra work to set up a define_external. That's why I posted a link to an example of doing so. And that's why I asked whether people would be interested in a Maple Library package that contains pre-made Maple (Library, interpreted) language level entry points to the routines.

This topic has only come up two or three times in this forum's history (another one being a request for LDL decomposition, many years ago I think, to which I also answered with the code for a wrapperless define_external to a CLAPACK function.) It's often hard to tell, when people don't ask for something much, whether they don't want it or whether they don't realize it might be asked.

ps. Are any of the CARE/DARE/LyaponovSolve/SylvesterSolve routines possible alternatives to the kind of DS system you need solved?

acer

@TomM You wrote, "However, the whole package is not bundled with Maple, in the sense that the code (dlls) comes with the Maple program distribution." Are you claiming that the dll is shipped with Maple, or not? A dynamic shared CLAPACK library object is, in fact, included in the Maple installation on at least 32 and 64 bit MS-Windows and Linux.

The real and complex double precision funtions (Dxxxxx and Zxxxxx, but not single precision Cxxxxx and Sxxxxx) of CLAPACK 3.0 should be in the clapack.dll that is indeed shipped with Maple. If there are some functions in the dll which are not exported, then that would be something to fix.

There isn't a dynamic BLAS library object shipped with Maple to which you can dynamically link, expect (ATLAS) on Linux.

I agree with you that there is some extra work to set up a define_external. That's why I posted a link to an example of doing so. And that's why I asked whether people would be interested in a Maple Library package that contains pre-made Maple (Library, interpreted) language level entry points to the routines.

This topic has only come up two or three times in this forum's history (another one being a request for LDL decomposition, many years ago I think, to which I also answered with the code for a wrapperless define_external to a CLAPACK function.) It's often hard to tell, when people don't ask for something much, whether they don't want it or whether they don't realize it might be asked.

ps. Are any of the CARE/DARE/LyaponovSolve/SylvesterSolve routines possible alternatives to the kind of DS system you need solved?

acer

A natural followup here is to try and make the resulting exported image file show as more like a Maple plot, with axes or text or labels, etc.

This can be done by creating an empty plot with all the desired text and axes, exporting that to an image, reading the image back in, and then using that to mask the density image (which pretty much amounts to overlaying them, if the text is black).

An exported empty plot has a blank border, as the axes don't lie right at the image boundary. This can be accomodated by either fitting the density image to the axes-bounded portion of the empty plot's image, or vice-versa by chopping off the blank border or a slightly larger empty plot. I take the second approach, below. (I suspect the first would be a better choice, if one were aiming to produce a good, fully fledged plot exporter or plot driver replacement.) The boder size also varies with the choice of Maple interface, eg, commandline or Standard GUI.

This calls `ArgPlot` as coded in the Post. This is rough, top-level stuff, but ideally it could be made into a more general procedure.

F:=z->`if`(z=0,Pi,exp(1/z)):
hsize,vsize:=300,300:

st,str,ba,bu:=time(),time[real](),kernelopts(bytesalloc),kernelopts(bytesused):
  Q:=ArgPlot(F,-0.25,0.25,vsize,hsize):
time()-st,time[real]()-str,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

st,str,ba,bu:=time(),time[real](),kernelopts(bytesalloc),kernelopts(bytesused):
_Env_plot_warnings:=false:
tempfile:=cat(kernelopts('homedir'),"/temp.jpg"):
# The Standard GUI's plot driver pads the plot with whitespace.
# So we increment the target size for export, and then later chop.
# The padding seems to be absolute in size, rather than relative
# to the dimensions, at least for the GUI's jpeg plot device. Certainly the
# padding can be much greater for jpeg export using the commandline interface's
# jpeg plot driver. Adjust accordingly.
plotsetup('jpeg','plotoutput'=tempfile,'plotoptions'=cat("height=",vsize+6+4,
                                                         ",width=",hsize+5+5));
plot(I,-0.25..0.25,'view'=[-0.25..0.25,-0.25..0.25]);
plotsetup('default'):
try fflush(tempfile); catch: end try:
fclose(tempfile);
# Why doesn't the GUI plot driver flush before returning control?!
# The timing is affected by a delay (YMMV, on NFS or a network drive)
# to wait for the image file to be written. If the `Read` below fails
# and claims the file does not exist, the delay may need increasing.
Threads:-Sleep(1): # 1 second
_Env_plot_warnings:='_Env_plot_warnings':
T:=ImageTools:-Read(tempfile,'format'='JPEG'):
correctedT:=rtable_redim(T[1+6..-1-4,1+5..-1-5,1..-1],1..vsize,1..hsize,1..3):
# The mask could likely be done more efficiently in-place on Q ("by hand"),
# with saving on allocation and thus subsequent collection.
newQ:=ImageTools:-Mask(Q,correctedT):
ImageTools:-Write(cat(kernelopts('homedir'),"/arg_expinv.jpg"),newQ,'quality'=95):
time()-st,time[real]()-str,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

ImageTools:-View(newQ);

 

This can be compared with less attactive results using `densityplot`. (I understand that densityplot's `grid` option is not the same as the width/height of an image. But increasing the `grid` values, for densityplot or complexplot3d, can be prohibitive in terms of time & memory resources and the performance of the plot driver/exporter and of the GUI itself. These aspects were the original motivation. The usual plotting commands cannot match the image creation method here, especially at higher or publication quality resolutions.)

restart:

F:=z->`if`(z=0,Pi,exp(1/z)):
hsize,vsize:=300,300:
plots:-densityplot(-argument(exp(1/(a+b*I))),a=-0.25..0.25,b=-0.25..0.25,
                   colorstyle=HUE, style=patchnogrid,scaletorange=-Pi..Pi,
                   grid=[hsize,vsize],axes=normal);

Exporting using context-menus (right-click, takes a long time...) produces this,

And using the GUI's plot driver with the plotsetup command produces this,

plotsetup('jpeg','plotoutput'="driv_expinv.jpg",'plotoptions'=cat("height=",300,
                                                                  ",width=",300));

acer

As far as the computational part goes, there seems to be some bottleneck in mapping log[2] over the complex-valued Matrix D2. Ideally, that subtask should be solvable with a very fast implementation using the Compiler, but I've put in comments about some weirdness with the Compiled procedure's behaviour there. (I'll submit an bug report.)

For producing higher resolution graphics than what densityplot (with colorstyle=HUE) provides, you might have a look at this recent Post. At higher resolutions, building the image directly can be much faster and more memory efficient than using densityplot. (The GUI itself uses huge amount of unreported memory for larger grids.)

The code below has some tweaks, with original lines commented out. (Uncommenting and running side by side might give some insight into why some ways are faster.) My machine runs the code below, with n=2^8, in under 7 seconds, and the originally Posted code in about 18 seconds. But I believe that a well-behaved Compiled complex log[2] map over a complex[8] Matrix would allow it to be done in around 1 second.

If you want the computational portions much faster, you could consider going with Compiled procedures to do all the computation, and give up the `~` mapping and evalhf. That also ensures sticking with hardward datatype Matrices (which didn't happen with the original code) and helps keep the memory use down.

restart:
with(DiscreteTransforms): with(plots):

n:=2^8:

Start:=time():

M := Matrix(n, n, 1, datatype = float):
Ii := Vector(n, proc (i) options operator, arrow; i end proc, datatype = float):
x := Vector(n, proc (i) options operator, arrow; Ii[i]-(1/2)*n end proc):
y := Vector(n, proc (i) options operator, arrow; (1/2)*n-Ii[i] end proc):
X, Y := Matrix(n, n, proc (i, j) options operator, arrow; x[j] end proc),
        Matrix(n, n, proc (i, j) options operator, arrow; y[i] end proc):
R := 10:

st:=time():
A := map[evalhf](`<=`,`~`[`+`](`~`[`^`](X, 2), `~`[`^`](Y, 2)), R^2):
time()-st;

#st:=time():
#AA := `~`[`@`(evalb, `<=`)](`~`[`+`](`~`[`^`](X, 2), `~`[`^`](Y, 2)), R^2):
#time()-st;

st:=time():
M := evalhf(zip(`and`,(A, M))):
time()-st;

#st:=time():
#MM:=`~`[evalhf](`~`[`and`](A, M)):
#time()-st;
#LinearAlgebra:-Norm(M-MM);

#plots:-listdensityplot(M, smooth = true, axes = none);

D1 := FourierTransform(M, normalization = none):
D2 := ArrayTools[CircularShift](D1, (1.0/2)*n, (1.0/2)*n):

st:=time():
Q1:=Matrix(map[evalhf](t->Re(abs(t)),D2),datatype=float[8]):
time()-st;

#st:=time():
#Q11:=abs(D2);
#time()-st;
#LinearAlgebra:-Norm(Q1-Q11);

st:=time():
plots:-listdensityplot(Matrix(map[evalhf](t->Re(abs(t)),D2),datatype=float[8]),
                       smooth = true, axes = none, colorstyle = HUE,
                       style=patchnogrid);
time()-st;

st:=time():
D3:=`~`[subs(LN2=evalhf(ln(2)),t->ln(t)/LN2)](D2):
time()-st;

#st:=time():
#D33 := `~`[log[2]](D2):
#time()-st;
#LinearAlgebra:-Norm(D3-D33);

# There is a very weird behaviour, where both evalhf and a Compiled version of
# this makeD3 routine can (for n=8 say) take either less than a second, or 25 seconds,
# intermittently. The very same calls, without re-Compilation, but vastly different timing.
# I've commented the attempts out, and am using the uncommented log[2] method above.
makeD3:=proc(n::integer,a::Matrix(datatype=complex[8]),b::Matrix(datatype=complex[8])) local i::integer, j::integer, LN2::float; LN2:=1.0/ln(2.0); for i from 1 to n do for j from 1 to n do #b[i,j]:=log[2](a[i,j]); b[i,j]:=ln(a[i,j])*LN2; end do; end do; NULL; end proc: #st:=time(): #cmakeD3:=Compiler:-Compile(makeD3,optimize=true): #time()-st; #st:=time(): #D3:=Matrix(n,n,datatype=complex[8]): #evalhf(makeD3(n,D2,D3)): #time()-st; #D3:=Matrix(n,n,datatype=complex[8]): #st:=time(): #cmakeD3(n,D2,D3); #time()-st; st:=time(): plots:-listdensityplot(abs(D3), smooth = true, axes = none, colorstyle = HUE, style=patchnogrid); time()-st; time()-Start;

[note: every time I edit this code it screws up the less-than synbols in the pre heml tags and misinterprets, removing content.]

aperture_density.mw

acer

As far as the computational part goes, there seems to be some bottleneck in mapping log[2] over the complex-valued Matrix D2. Ideally, that subtask should be solvable with a very fast implementation using the Compiler, but I've put in comments about some weirdness with the Compiled procedure's behaviour there. (I'll submit an bug report.)

For producing higher resolution graphics than what densityplot (with colorstyle=HUE) provides, you might have a look at this recent Post. At higher resolutions, building the image directly can be much faster and more memory efficient than using densityplot. (The GUI itself uses huge amount of unreported memory for larger grids.)

The code below has some tweaks, with original lines commented out. (Uncommenting and running side by side might give some insight into why some ways are faster.) My machine runs the code below, with n=2^8, in under 7 seconds, and the originally Posted code in about 18 seconds. But I believe that a well-behaved Compiled complex log[2] map over a complex[8] Matrix would allow it to be done in around 1 second.

If you want the computational portions much faster, you could consider going with Compiled procedures to do all the computation, and give up the `~` mapping and evalhf. That also ensures sticking with hardward datatype Matrices (which didn't happen with the original code) and helps keep the memory use down.

restart:
with(DiscreteTransforms): with(plots):

n:=2^8:

Start:=time():

M := Matrix(n, n, 1, datatype = float):
Ii := Vector(n, proc (i) options operator, arrow; i end proc, datatype = float):
x := Vector(n, proc (i) options operator, arrow; Ii[i]-(1/2)*n end proc):
y := Vector(n, proc (i) options operator, arrow; (1/2)*n-Ii[i] end proc):
X, Y := Matrix(n, n, proc (i, j) options operator, arrow; x[j] end proc),
        Matrix(n, n, proc (i, j) options operator, arrow; y[i] end proc):
R := 10:

st:=time():
A := map[evalhf](`<=`,`~`[`+`](`~`[`^`](X, 2), `~`[`^`](Y, 2)), R^2):
time()-st;

#st:=time():
#AA := `~`[`@`(evalb, `<=`)](`~`[`+`](`~`[`^`](X, 2), `~`[`^`](Y, 2)), R^2):
#time()-st;

st:=time():
M := evalhf(zip(`and`,(A, M))):
time()-st;

#st:=time():
#MM:=`~`[evalhf](`~`[`and`](A, M)):
#time()-st;
#LinearAlgebra:-Norm(M-MM);

#plots:-listdensityplot(M, smooth = true, axes = none);

D1 := FourierTransform(M, normalization = none):
D2 := ArrayTools[CircularShift](D1, (1.0/2)*n, (1.0/2)*n):

st:=time():
Q1:=Matrix(map[evalhf](t->Re(abs(t)),D2),datatype=float[8]):
time()-st;

#st:=time():
#Q11:=abs(D2);
#time()-st;
#LinearAlgebra:-Norm(Q1-Q11);

st:=time():
plots:-listdensityplot(Matrix(map[evalhf](t->Re(abs(t)),D2),datatype=float[8]),
                       smooth = true, axes = none, colorstyle = HUE,
                       style=patchnogrid);
time()-st;

st:=time():
D3:=`~`[subs(LN2=evalhf(ln(2)),t->ln(t)/LN2)](D2):
time()-st;

#st:=time():
#D33 := `~`[log[2]](D2):
#time()-st;
#LinearAlgebra:-Norm(D3-D33);

# There is a very weird behaviour, where both evalhf and a Compiled version of
# this makeD3 routine can (for n=8 say) take either less than a second, or 25 seconds,
# intermittently. The very same calls, without re-Compilation, but vastly different timing.
# I've commented the attempts out, and am using the uncommented log[2] method above.
makeD3:=proc(n::integer,a::Matrix(datatype=complex[8]),b::Matrix(datatype=complex[8])) local i::integer, j::integer, LN2::float; LN2:=1.0/ln(2.0); for i from 1 to n do for j from 1 to n do #b[i,j]:=log[2](a[i,j]); b[i,j]:=ln(a[i,j])*LN2; end do; end do; NULL; end proc: #st:=time(): #cmakeD3:=Compiler:-Compile(makeD3,optimize=true): #time()-st; #st:=time(): #D3:=Matrix(n,n,datatype=complex[8]): #evalhf(makeD3(n,D2,D3)): #time()-st; #D3:=Matrix(n,n,datatype=complex[8]): #st:=time(): #cmakeD3(n,D2,D3); #time()-st; st:=time(): plots:-listdensityplot(abs(D3), smooth = true, axes = none, colorstyle = HUE, style=patchnogrid); time()-st; time()-Start;

[note: every time I edit this code it screws up the less-than synbols in the pre heml tags and misinterprets, removing content.]

aperture_density.mw

acer

@newlam All I see is base64 encoded source in your post, not images.

You can no longer paste 2D Math from Maple straight into Mapleprimes, since the site was revised. But mere images of your code can be awkward for us, since we don't want to have to type out anything.

You can use the green up-arrow in the editor here, to upload your worksheet. Or you can paste the code as 1D text Maple input.

Forcing the tickmarks seems to be key to this kind of thing, in current Maple. And that can be done for the x-axis or y-axis separately, or for both.

Flipping the plot in the x-direction can often be accomplished by a simple change of variables in the expression (as in the link Kamel has given). But it can also be done for either the x- or the y-direction using plottools, eg here.

acer

Forcing the tickmarks seems to be key to this kind of thing, in current Maple. And that can be done for the x-axis or y-axis separately, or for both.

Flipping the plot in the x-direction can often be accomplished by a simple change of variables in the expression (as in the link Kamel has given). But it can also be done for either the x- or the y-direction using plottools, eg here.

acer

First 433 434 435 436 437 438 439 Last Page 435 of 592