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

Are you trying to say that you want to avoid the overhead of three calls to Transpose (or any kind of copying)? I mean, are you trying to avoid any overhead incurred by, say,
    Transpose(LinearSolve(A^%T,b^%T),inplace)

How large is your Matrix? Are the entries all floats? Is b always a Vector? Does the Matrix have an indexing-function (ie, shape)?

Are you looking for a convenient syntax, or is high efficiency a primary concern?

In the floating-point case, are you concerned with possibly less accuracy from b.A^(-1)?

Here is a simpler example with only a single plot rather than an animation.

In the Classic interface of Maple 2019.0 (32bit Windows) the PLOT3D structure that contains two LIGHT substructures will render using both light sources together.

In the Standard GUI the renderer only uses one of the light sources. As far as I know this has always been the case.

When the plotting commands were retro-fitted to use the new keyword handling (?paramprocessing) any arguments of the form light=[...] got the new handling so that only the last instance gets utilized. This change came some time after Maple 7. Prior to that multiple instances of the light=[...] option were all utilized and the resulting PLOT3D structure contained the corresponding multiple LIGHT substructures (but only Classic would utilize them in rendering).

In summary:
1) The Standard GUI renders using only one of any LIGHT substructures found, even when there are several.
2) The Classic GUI (32bit Maple 2019 on MS-Windows) can still render using multiple LIGHT substructures, if they are present.
3) The parameter-processing of the plotting commands now ignores all but the last light=[...] keyword option, so there no longer is any convenient way to construct a 3D plot that contains multiple LIGHT substructures.

Here's an example that I did in Maple 2019.0, 32bit Classic on Windows 7. I have replaced the images, so that it renders here as it does in Classic. (In the Standard GUI the first plot, which contains both LIGHT substructures, renders the same as the third plot. The default orientation may differ from Classic, but that's not the point.)

restart;
with(plots): with(plottools):
P := display(sphere([0,0,0],1,color=gray,style=surface)):
L1:=LIGHT(90, 40, 1, 0.5, 0.8):
L2:=LIGHT(90, -80, 1, 0.2, 0.1):

op(0,P)(op(P),L1,L2);


op(0,P)(op(P),L1);


op(0,P)(op(P),L2);


 

Download multiplelights.mws

@sand15 I recommend getting Maple 2019.

It comes bundled with LLVM.

For the very best experience in terms of stability and performance I'd also recommend it (or any modern version) on Linux.

@mmcdara For fun, here using Threads:-Seq and a few other storage twists.

Here, 602.30ms for the Box-Muller versus 466.70ms for the default Sample method and 2.31s for Sample's envelope method, in real-time, at N=10^7.

CodeTools:-Usage(Statistics:-Sample(Normal(0,1), 10^7, method=envelope), iterations=10):
memory used=95.36MiB, alloc change=0.52GiB, cpu time=2.33s, real time=2.31s, gc time=21.52ms

CodeTools:-Usage(BoxMuller_ac3(10^7), iterations=10):
memory used=79.96MiB, alloc change=0.78GiB, cpu time=938.50ms, real time=602.30ms, gc time=0ns

CodeTools:-Usage(Statistics:-Sample(Normal(0,1), 10^7), iterations=10):
memory used=76.31MiB, alloc change=0.75GiB, cpu time=466.30ms, real time=466.70ms, gc time=0ns

I was disappointed with the speedup under Threads:-Seq, though I did not try to find the optimal splitting length, given that there is overhead that can swamp if too many threads are allowed for tasks that are too small.

Using Threads:-Task and binary splitting seems a bookkeeping headache, although perhaps a spot in each subvector could store its length and the number of relevant entries.

Also, as written these procedures could actually return more than N values, although it would not be much effort to restrict this efficiently to at most N.

[edit] Apparently I forgot to attach the revision. I'll do that later today.

[edit] I also realize that another 5-10% of the time (and some memory overhead) can be shaved off by forming result R with exactly N entries, and quitting the loop early once N results are attained, etc. I will revise this evening.

@Carl Love I think that you're right, in the sense that it's not clear whether the OP's example is input or the result of some earlier computations. Since the OP asked just last week about substituting into the result of a series call, this Question is less clear than it could be.

@mmcdara Ok, so you are getting 64ms for N=10^5 with _ac2 under evalhf mode. That's a good start.

On 64bit platforms recent versions of Maple are supposed to use the LLVM compiler by default (unless the inmem=false options is passed to Compile, in which case it tries to use gcc from the OS). The LLVM compiler is bundled with the installation of Maple itself.

If your Maple version is "older" you might have to install gcc on your OSX (or Linux) box. It's free.

I don't recall which was the first Maple version to ship with LLVM bundled on the 64bit installs.

If you cannot get the Compiler to work in any version, including 2019.2 say, then there must be some internal problem that prevents it from working. It'd be nice to get this resolved for you.

Mac Dude showed you convert(...,polynom) in his Answer to your series Question from one week ago.

@dharr 

Yes, that is an alternate, convenient syntax for obtaining a procedure T which (just like in my second worksheet) returns unevaluated for nonnumeric argument c.

The overall time to compute can be reduced by about a factor of 2 by giving it option remember (just as in my revision of my second worksheet). That is,

   T := unapply('fsolve'(m(t, c) = -1, t = 0 .. 1/2*Pi), c, numeric, proc_options = [remember])

This is a crude revision, without any effort to handle the excluded entries better.

(I'm in a hurry just now. Later, I'll look at it more carefully. Please don't consider this as mu "submission". It's just to show some simple ideas, keeping the natural style somewhat... Quite a bit more optimization is possible.

But there is also an X1,X2 aspect that I'm not content with.)

BoxMuller_ac1 := proc(N)
  local V, S, X1, X2:
  V  := Statistics:-Sample(Uniform(-1, 1), [ceil(2*N/3), 2]):
  S  := V[..,1]^~2 +~ V[..,2]^~2;
  map[evalhf,':-inplace'](proc(t)
                           if t>1.0 then return Float(infinity); end if;
                           (-2*log(t)/t)^(1/2);
                         end proc,S);
  X1 := select[':-flatten'](type, S *~  V[..,1], numeric);
  X2 := select[':-flatten'](type, S *~  V[..,2], numeric);
  return < X1 , X2>:
end proc:

Tell us the URL for the website of this course, so that we can download the complete details, including the definitions of those procedures.

@mmcdara 

A considerable speedup can be had from using just evalhf (as I showed). True, the Compiler brings more (as I also showed), but evalhf can be quite handy.

The Compiler is still bundled and working in my Maple versions 2018 and 2019.

The is command is not builtin, and is not fast. Even use of a custom anonymous operator as the first argument to select or Select means you've moved away from builtins. Pure builtin use would more like, say,

     select(`>`,S,q)

When you do both Heaviside and subtraction as separate elementwise operations you produce a garbage Vector of length N for each operation. Overhead costs add up.

Apart from your concerns with the behavior of Sample (with its default method), I am surprised by the performance of add on a float[8] Vector.

In the examples below the only difference is the approach to adding up the numbers of outliers in the generated samples. The results concur, but the performance varies.

I would have expected that the first of these would perform comparably to the second.

(I include the others for your interest.)

restart;

with(Statistics):

X := RandomVariable(Normal(0, 1)):

N := 10^6: R := 20: Q := [$0..5]: nQ := numelems(Q):

M  := Matrix(R, nQ, 0):

str:=time[real]():
for r from 1 to R do
  S := Sample(X, N):
  for q in Q do
    M[r, q+1] := add(Heaviside~(S -~ q)):
  end do:
end do:
time[real]()-str;

34.479

restart;

with(Statistics):

X := RandomVariable(Normal(0, 1)):

N := 10^6: R := 20: Q := [$0..5]: nQ := numelems(Q):

M  := Matrix(R, nQ, 0):

str:=time[real]():
for r from 1 to R do
  S := Sample(X, N):
  for q in Q do
    tt := Heaviside~(S -~ q):
    M[r, q+1] := evalhf(add(tt[i],i=1..N)):
  end do:
end do:
time[real]()-str;

6.561

restart;

with(Statistics):

X := RandomVariable(Normal(0, 1)):

N := 10^6: R := 20: Q := [$0..5]: nQ := numelems(Q):

M  := Matrix(R, nQ, 0):

K1:=proc(ss,q,N) local i; add(ss[i]>q,i=1..N); end proc:

str:=time[real]():
for r from 1 to R do
  S := Sample(X, N):
  for q in Q do
    M[r, q+1] := evalhf(K1(S,q,N));
  end do:
end do:
time[real]()-str;

5.800

restart;

with(Statistics):

X := RandomVariable(Normal(0, 1)):

N := 10^6: R := 20: Q := [$0..5]: nQ := numelems(Q):

M  := Matrix(R, nQ, 0):

K2:=Compiler:-Compile(proc(ss::Vector(datatype=float[8]),q::float,N::integer)
                        local i; add(`if`(ss[i]>q,1,0),i=1..N);
                      end proc):

str:=time[real]():
for r from 1 to R do
  S := Sample(X, N):
  for q in Q do
    M[r, q+1] := K2(S,q,N);
  end do:
end do:
time[real]()-str;

1.655

 

Download mmc_comp.mw

That was Maple 2019.1 on 64bit Linux. I also tried Maple 2017.2, in which the relative performance was similar.

@Carl Love That's also what I was trying to drive at.

I'd rather see the original problem then some concocted (series?) approximation. And the use of limit as y->infinity on the pade approximation is also odd. The whole thing looks shaky.

@Puneeth Anyone can observe your final call to solve(...,[a,b]) and figure out that you're trying to find a and b. 

What does your computation represent? What do the expressions (which you pass to `pade`) represent? What does n=11 represent?

One mistake is that you set Digits too low. Another is that you use `sum` in some places where `add` would be better.

Another is that (in the statement that assigns to F[k+3] in the loop) you have an instance of F(k+2) which is likely a typo for F[k+2].

Another is that you're missing the closing brace } in the call to solve.

But a bigger mistake is in not telling us properly what you're trying to accomplish.

First 197 198 199 200 201 202 203 Last Page 199 of 592