Why do you need to call `solve` more than once here?

restart:

tau:=proc(jj)
   if not type(jj,numeric) then
      return 'procname'(args); end if;
   cat(tau,jj); end proc:
eta:=proc(cc)
   if not type(cc,numeric) then
      return 'procname'(args); end if;
   cat(eta,cc); end proc:
mix:=proc(jj,cc)
   if not ( type(jj,numeric) and type(cc,numeric) ) then
      return 'procname'(args); end if;
   cat(mix,jj,cc); end proc:

basesol:=solve({log(pC[j,c]/(1-pC[j,c]))=mu+tau(j)+eta(c)+mix(j,c)},pC[j,c])[1]:
solform:=eval(pC[j,c],basesol):

K:=50;C:=20;
cat(p,'C',1..20):=seq([seq(solform,j=2..K)],c=1..20):

acer

I believe that writing it from scratch in Maple would be an easier way to get an efficient implementation than would be attempting as you have here to translate (close to line-by-line) an implementation from Matlab or some other language.

I particular, the many calls to SearchArray (to replicate `find`, say), the many elementwise `~` operations, all producing copies of Arrays, is slow. Not using datatype=float[8] throughout does not help.

I suggest doing initial edge detection -- using convolution say. Then threshold a bit, perhaps. The write a proc to populate a float[8] accumulator Array, inplace, using the edge info Array. Make it evalhfable & Compilable, Then right another proc which analyzes the accumulator results (to "de-Hough" or threshold, etc). The at the very end, a proc which draws the deduced circles/shapes on an image Array. This should be able to run in a few seconds. And one of the expensive bits -- the accumulator population -- can likely be parallelized by Task splitting across all the deteced edge points.

I'm excited to try this approach. But not for a few days at the very least...

acer

Maybe you can just scale the image down, and then do the 3D matrixplot with a few modifications?

restart;
with(ImageTools):
with(plots):

a:=Read(cat(kernelopts(homedir),"/lighthouse1.jpg")):

c:=Scale(ToGrayscale(a),1..300):
rtable_options(c,'subtype'=Matrix); # no need for double copying

matrixplot(map[evalhf](`*`,c,200),style=point,shading=zhue,orientation=[0,0],
           symbolsize=5,scaling=constrained);

The image is simply scaled down, to get better GUI responsiveness. The scaling=constrained option better respects the width/height ratio. And by muliplying the z-values (above, by 200) you can get more depth visible when rotated in 3D. And by adjusting the symbolsize you may gain some effects (eg. on the bleedthrough). So you can vary and experiment with those three adjustments.

I just tried it in Maple 12 as well. Such amazing performance for the 3d plot rotation! (I'm using 64bit Linux and an Nvidia card, with hardware accel. purportedly enabled. But it is much slower in Maple 16/17.)

The color intensity is much nicer too, on the 3D plot's zhue shading. But that can be obtained in 16/17 with the lightmodel=none option.

acer

The Standard GUI is relatively slow to render and rotate 3D density-style plots, if the number of grid points is moderately large.

I can't quite tell whether making a rotatable 3D plot is your actual goal, or not...  If it is, then you'll have to reduce the image somehow, sure.

Another possiblity for a 3D plot might be to make use `surfdata` and subsop the image Array into into COLOR substructure. But you'd likely still need to reduce the image somehow. This would give a smoother surface than would `matrixplot`, but maybe you want the granularity of the latter?

BTW, the values of the grayscale conversion are just one way amongst several possible for getting the "intensity". Another way is to grab the 3rd ("V"=intensity) layer of the HSV conversion. And since you seem to want zhue shading, then one way to get that is to replace the 1st ("H"=hue) layer with the 3rd layer after doing a RGBtoHSV conversion.

This attachment reads an image file, does the replacement of "H" layer by the "V" layer (force-spread to 0..1), and writes out the result after max'ing out both the "S"=saturation and "V"=intensity layers. Then it sticks both on a pair of Labels. Sure, you can't rotate it in 3D, but it works pretty fast and the GUI still behaves well afterwards. This approach is more of a replacement for 2D `densityplot` or `listdensityplot` (for which the Standard GUI responds quite a bit worse than it does for even the 3D plot). You may, or may not, want to force-spread the hue to match 0..1 end-points.

[edit: You could get fancy with the Labels and use SetProperty to adjust the Label componet's pixelHeight and pixelWidth to be in the same propertion as the image Arrays first two dimensions, and even scale down the image Array to match the number of pixels displayed on the Label, to minimize file i/o). I didn't do all that.]

VasH.mw

acer

I would have expected to see something within 0.6 ulps of 7.100000001e8 and 7.0999999820e8 seems too far off for the default working precision setting Digits=10.

This is already better,

restart: 
Digits:=11:
2^29.403243784;

                                       8
                        7.1000000014 10 

After all, is this not an "atomic" operation?

Here's something fun, for all the lovers of automatic simplification. Compare the results of these various computations (the restarts are necessary). The middle one is the one to note.

restart:
a:=29.403243784:
for i from 10 to 50 do
  Digits:=i:
  2^a;
end do;

restart:
for i from 10 to 50 do
  Digits:=i:
  2^29.403243784;
end do;

restart:
p:=proc() option trace; local i;
for i from 10 to 50 do
  Digits:=i:
  2^29.403243784;
end do;
NULL;
end proc;
p();

And if you really love automatic simplification then you might like these too:

restart:

4^(1/2);

                              (1/2)
                             4     

'(4^(1/2))^29.403243784';

                                      8
                        7.099999982 10 

2^(1/2), type(2^(1/2),numeric), type(Pi,numeric);

                       (1/2)              
                      2     , false, false

'(2^(1/2))^29.403243784', 'Pi^29.403243784';

                                 29.403243784
                  26645.82515, Pi            

Since exact 2^(1/2) is not rational it might seem a bit weird for that example to auto-simplify.

acer

There is a line,

alpha=evalf(alpha);

Making it into an assignment appears to work. At the beginning, alpha is aliased to LambertW(exp(-1)) which is exact.

But perhaps exact LambertW(exp(-1)) is a problem during the numeric dsolving...  Does someone see a better way, which allows instances of alpha to remain with its symbolic appearance throughout? I mean, something more graceful than either assigning alpha as a float or substitution of alpha=evalf(alpha) into each or the problematic DE system arguments to DEplot?

acer

Your code has this incorrect syntax,

  with*Optimization;

which should probably be

  with(Optimization);

Also, if that's the package you are trying to use then the command within it is `Minimize` not `minimize`. Note that Maple is case-sensitive.

If you fix those two items, and still have problems with it, then please just add a followup comment in this current thread.

acer

It might be tempting to believe that the answer could be (2*n*2^a*Pi^(a+1)+2*a+2)/((a+1)*n).

And a hint might be seen by the crude substitution of values,

restart:
seq( [int(abs(sin(n*x)-x^i),x=0..2*Pi),
      normal(2/(i+1)*(2^i*Pi^(i+1)*n+i+1)/n)],
     i=2..10 ) assuming n::posint;

But there is a difference of 2/n whioh I don't understand in the following,

restart:

unprotect(series): oldseries:=eval(series):
unprotect(limit): oldlimit:=eval(limit):
series:=MultiSeries:-series:
limit:=MultiSeries:-limit:

cand:=int(abs(sin(n*x)-x^a),x=0..2*Pi) assuming a>1, n::posint;

    1     /   (a + 1)   (a + 1)   /     /      /             a\ 
--------- \n 2        Pi        - \limit\signum\-sin(n x) + x / 
n (a + 1)                                                       

                     /             a\  (a + 1)  
  cos(n x) a + signum\-sin(n x) + x / x        n

           /             a\                       \\        \
   + signum\-sin(n x) + x / cos(n x), x = 0, right// + a + 1/

series:=eval(oldseries):
limit:=eval(oldlimit):

simplify(cand) assuming a>1, n::posint;

                        (a + 1)   (a + 1)
                       2        Pi       
                       ------------------
                             a + 1   
    
simplify( (2/(a+1)*(2^a*Pi^(a+1)*n+a+1)/n) - % );

                               2
                               -
                               n

I suppose that it might not necessarily be valid to apply (and simplify under) an assumption like n::posint to `cand` above while it still contains a pending unevaluated limit (as x->0+).

acer

Explore(plot(A*sin(w*x+c),x=-2*Pi..2*Pi,view=-10..10),
        parameters=[A=1.0..10,w=1.0..20,c=Pi/6..11*Pi/6]);

You could also do,

plots:-interactive(A*sin(w*x+c));

and use the uppermost drop-menu in the PlotBuilder dialogue which pops up, and select interactive plots with one/two parameters.

You can also use `plots:-interactiveparams` and get a Maplet popup view with sliders.

And some of these can also work by calling `plots:-animate` as their plotting command.

See the help-pages for all these commands, for more examples.

acer

There are several things that would make life easier here.

Use `*` on your keyboard for multiplication, even if you are in 2D Math input mode. Using that cross symbol (from a palette, presumably) gets something that evaluates ok, sure. But it is ugly, and conversion to 1D Notation goes pretty wrong as Patrick mentioned.

Your lowercase letter L displays almost exactly like the numeral one, which makes it very confusing visually and leads to visual ambiguity (and perhaps even mistakes). For example, there is no lowercase L appearing in the integrand of your second expression. Is that right? Hard to tell what you intended. Oh, wait, that upper limit of integration is... yes.. a lowercase L and not the number 1. It's confusing, since on my Windows 7 Maple 17 those  prettyprint almost identically for your input expression. So I will use capital L from here on.

Next, for unknown `d` and `L` your integral will not compute exactly (well, it might under assumptions, but let's keep it simple for now). This is the second time in 24 hours that this topic has arisen: plotting an integral using the active `int` command which fails to evaluate exactly can be slow. It can be slow because -- after the integration of the argument involving active `int` fails  -- the plotting routine receives and unevaluated `int` call as its argument. But then, for every substituion of numeric values of the plotting variables, this integral gets attempted over and over. Now, you might get lucky if your integration range also happens to be numeric (since a float range will make `int` do numeric quadrature), but really the whole topic is too confusing for the new user at this point.

... so it's often much simpler and faster to use inert uppercase `Int` for integral expressions that are to be plotted.

Another possibility, which I do not show below, is to try and get the exact value of the integral under assumptions on `d` and `L` which match the plotting ranges you've supplied.

Now, we've seen a hint that you may be relying on the palettes for entering 2D Math expressions (as an unfortunate consequence of over-enthusiam in the help portal for Clickable Math, perhaps). And the Expression palette only provides the active (dark black integral symbol) `int`. Well, you can still get a nicely typeset inert integral `Int` by using command-completion. Just start typing Int as 2D Math input, and then do  command-completion (Ctrl-spacebar on MS-Windows.) The second entry in the menu the pops up will be the template for an inert definite integral.

Notice that inert `Int` gets typset in a lighter shade of gray for the integration symbol.

There is another possible reason why you might want an inert integral. Your plotting and integration variables take on values in the region near 10^(-6), and the integrand takes 4th powers in the denominator and scales  heavily by numeric factors in the numerator. This might be bad scaling. Which might be why you encountered some numeric difficulties, ...but then I didn't really test for that. Anyway, if you want to then you can rescale either `d` or `L` or both. It's quite possible that your original is better scaled, and better left as is. I really didn't check.

Here's a worksheet that might help.

If you decide to rescale then you may wish to use axis labels which show the scaling factors. At the end of the sheet are a few different ways to get such products to be typset. None is the (perfect, IMHO) visual form like `d 10^-6`, however.

Download dropmodif.mw

acer

Does this make something like what you're after? (It requires m>=2, n>=2 which you could put in as argument checking if you have the need.)

G:=proc(m,n)
   local minmn, maxmn;
   minmn:=min(m,n); maxmn:=max(m,n);
   LinearAlgebra:-BandMatrix([[1$`if`(n>=m,minmn-2,n-1),2],
                              [3$minmn],
                              [2,1$`if`(n>m,m-1,minmn-2)]],1,m,n);
   end proc:

G(5,4);
G(10,7);
G(2,2);
G(6,6);
G(4,5);
G(7,10);

If it's not quite what you want, in the case of nonsquare n<>m, then you could try adjusting the code after looking over the BandMatrix help page. Note that this makes a `band[1,1]` storage Matrix. You can change that, or the datatype if you need it for speed efficiency, as options to that constructor command.

acer

This doesn't recover subsection 2.1, but it does seem to get back the earlier part of section Lek 2.

BETmodif.mw

acer

You seemed to be missing a closing bracket, and so I made a guess as to which variable you wanted as the dummy variable of integration. (Adjust, if I guessed wrongly.)

How about something like this? It forces purely numerical integration and restrains its requested accuracy because you might not need so many digits of accuracy, just for plotting.

restart:

Z := theta->sin(theta+2.3)+sin(2*theta+2*(-1.5))+sin(3*theta+3*(-0.94e-1)):
Zbar := (theta, phi)-> Z(theta+phi)-Z(theta):

H:=Int(signum(Zbar(theta, phi)+l)*abs(Zbar(theta, phi)+l)^(1/4),
       theta=0 .. 2*Pi, epsilon=1e-5, method=_d01ajc):

evalf(subs(phi=1/2,l=1/2,H));

plot3d(H, phi=0..1, l=0..1, grid=[20,20], axes=box);

acer

Two main things come to mind.

1) There are a few other ways to get similar results, some of which might already be acceptable (for some people).

For example, given your 20-by-20 Matrix `A`,

A[..,[seq(j, j=1..20, 2)]];

A[[seq(i, i=1..20, 10)],1..3];

2) Maple has some nice facilities for rtable indexing. But less is documented about being efficient, and yet it seems to me that quite often people are surprised when their easy-to-use, convenient code does not scale best as the problem sizes are ramped up high. If I had to make lots of such plots then I would be tempted to re-use structures rather than index into existing structures. For example I might create an m-by-2 float[8] Matrix `M` just once, and then ArrayTools:-Copy various rows or columns into `M` whenever I wanted to produce another such plot. It doesn't have to be done exactly that way, of course, and the methodology could reflect the situation. But for very large sizes I might be tempted by the idea of producing less work for the garbage collector, with each plot. Of course, such approaches might only be worth it if the interface (commandline, say) is also efficient in producing/exporting plots.

2) is also related somewhat to this (in which Joe's comment about ArrayTools:-Alias is worth noting.)

restart:
with(Statistics): 

P := Sample(RandomVariable(Normal(10, 10.)), 10^6):
F:=proc(x) `if`(ilog10(x)>-2,parse(sprintf("%.2f",x)),NULL); end proc:

X := RandomVariable(EmpiricalDistribution(map(F,convert(P,list)))):

plot(t->ProbabilityFunction(X,evalf[2](t)),-20..40);

plot(t->ProbabilityFunction(X,F(t)),-20..40);

acer