You should be using `eval` instead of `subs`, for the first way. In other words, eval(F,x=0) instead of subs(x=0,F) for example.

For your other way, you haven't acutally performed an "unapply" action. You've just written down an operator, manually. That will not work because the formal parameter of the operator (the x in x->) is not the same as the global x in the expression f. To get this way to work, you should do G:=unapply(lhs(f),x) instead.

acer

Here is a bit more.

2p81_final_modif2.mw

Perhaps you could make a backup copy.

The last version of updtsrc (suitable for getting MapleV source code file to run in Maple 6), can be obtained using the "HTTP" links from this page (which itself can be arrived at via a link from the ?updtsrc help-page).

But, isn't updtsrc a utility for upgrading plaintext source code only? I mean, it does language syntax converion, right? Does it also handle .mws? I thought not. If this is true, then would it work to open the old R4 .mws sheets in modern Maple's Classic GUI (or in MapleV R4's GUI, if you can), and then first to export (or copy & paste) the source code to plaintext .mpl files?

Where it gets sticky is with sheets older than R4, if I recall. There were some versions of updtsrc which could do R2->R3 conversion (and another for R3->R4 if I recall...). But the last one, linked to above, might not do those earlier conversions as well. I could be wrong about that. I have a hazy recollection of needing a multi-step process for making sheets for R2 or R3 become usable in Maple 7 or modern Maple. If anyone remembers this kind of detail, it'd probably be Thomas Richard.

If the version of updtsrc linked to from ?updtsrc doesn't fully work for the source code in your MapleV R4 sheets then let us know.

acer

First off, I think that one might need to set the Manipulator to "Click and Drag" so that that Action code is executed upon subsequent left-clicks, after the first animation. (I mention it, just in case anyone is having trouble running your example.)

If I am understanding you, then the issue is that the "current" frame of the old animation is displayed while the upcoming animation is being computed by the kernel and prepared for insertion by the GUI. I saw this more clearly if I used your example with frames=400. Then I ran it I stopped part way through, then slid the playbar back to the beginning or one side, then left-clicked to get the Action. What I saw was that the "last" frame of the previous animation was displayed while the upcoming animation was being prepared.

If I've understood rightly, then I should ask: what do you want displayed during that pause? The last shown frame of the previous animation, or a blank plot, or perhaps  the first frame of the upcoming animation? I think that any of these can be accomplished by adding in another line, before the insertion of the upcoming animation. Eg, uncommenting the one you want,

use DocumentTools, plots in

if cc = red then cc:= blue else cc:= red fi;

# quickly insert first frame of coming animation
#SetProperty(Plot0,value,plot(eval(sin(x+t),x=0),t=0..Pi,colour=cc),refresh=true);

# quickly insert blank frame
#SetProperty(Plot0,value,plot([[0,undefined]]),refresh=true);

# quickly insert something else
#...

SetProperty(Plot0,value,
            animate(plot,[sin(x+t),t=0..Pi,colour=cc],x=0..Pi,frames=400),
            refresh=true);
SetProperty(Plot0,play,true,refresh=true)

end use;

This is one of several reasons why I don't like using the traditional animation data structure inside Plot components. Another problem (for me) was that the animation runs asynchronously so any commands that come in the same code, but after the play command, get performed but then overridden by the animation which might still be playing asynchronously in its own thread. I was unable to find a reliable way to promptly apply an action that would occur only when the animation finished. The other main problem with traditional animations is, of course, their heavy use of memory resources in computing all frames up front. I believe that Maple is now fast enough that most "animations" can be done in an frame-by-frame computed-on-the-fly way, where all the controls would be done via library level & Component code. (Joe's comment about reusing a mere reference to an rtable inside a PLOT data structure may be extremely important here, for performance.)

acer

I'm not quite sure what you're asking. Are you unsure how to use the printf modifiers?

> s1:=0.5555552487:
> f(s1):=0.2542365647:

> printf( "s1 has value %.9e, and f(s1) has value %.10f\n", s1, f(s1) );

s1 has value 5.555552487e-01, and f(s1) has value 0.2542365647

See the printf help-page for more about the various modifiers (e, f, g, E, etc).

What do you mean by, "get it published"? If you want it printed to a file then use fprintf instead.

acer

Newton's method is what fsolve does for multivariate problems. You can see the source code for that interpreted Library routine in the usual ways, eg. by issuing the command,

    showstat(`fsolve/sysnewton`);

The source is a bit verbose, as it is adding guard digits, trying and be clever in the absence of a passed pair of stopping tolerances, handling both expression and operator form, etc.

> restart:

> f:=(a,b)->sin(a+b):
> g:=(a,b)->a-b:

> infolevel[fsolve]:=2:

> sol := fsolve({f,g});

fsolve/sysnewton: trying multivariate Newton iteration
fsolve/sysnewton:
guess vector [6.5647247568434774, -9.1035796650957585]
fsolve/sysnewton: norm of errors: 16.235204335176419
fsolve/sysnewton: new norm: 0.85321881380602701e-1
fsolve/sysnewton: iter = 1 |incr| = 15.668 new values _S01 = -1.6135 _S02 = -1.6135
fsolve/sysnewton: new norm: 0.20840812946897378e-3
fsolve/sysnewton: iter = 2 |incr| = 0.85634e-1 new values _S01 = -1.5707 _S02 = -1.5707
fsolve/sysnewton: new norm: 0.30173615373566167e-11
fsolve/sysnewton: iter = 3 |incr| = 0.20840e-3 new values _S01 = -1.5708 _S02 = -1.5708
fsolve/sysnewton: new norm: 0.38462643383279503e-16
fsolve/sysnewton: iter = 4 |incr| = 0.30174e-11 new values _S01 = -1.5708 _S02 = -1.5708
              sol := [-1.5707963267948966, -1.5707963267948966]

> f( op(sol) );

                                            -10
                              4.102067615 10   

> g( op(sol) );

                                     0.

kernelopts(version);
          Maple 13.02, X86 64 WINDOWS, Oct 12 2009, Build ID 436951

acer

You could have a look at this kind of thing (and also the various links at the bottom of that Comment, and others' Answers).

acer

The first argument, in Matrix form, will be like [-EV,2*Cov].

An example, with only simple bounds as constraints,

restart:
with(LinearAlgebra):
with(Optimization):

N:=600:
U:=RandomMatrix(N,outputoptions=[datatype=float[8],shape=triangular[upper]]):
H:=Matrix(U.U^%T,shape=symmetric):

EV:=RandomVector(N,outputoptions=[datatype=float[8]]):

bl:=Vector(N,fill=0,datatype=float[8]):
bu:=Vector(N,fill=1,datatype=float[8]):

sol1 := CodeTools:-Usage( QPSolve([-EV,2*H],[],[bl,bu]) ):
memory used=77.41MiB, alloc change=72.99MiB, cpu time=2.59s, real time=2.26s

W:=Vector(N,symbol=w):
sol2 := CodeTools:-Usage( QPSolve(W^%T.H.W-W^%T.EV,seq(W[i]=bl[i]..bu[i],i=1..N)) ):
memory used=243.44MiB, alloc change=50.99MiB, cpu time=5.05s, real time=5.06s

Norm(Vector(sol1[2])-Vector(eval(W,sol2[2])));
                               0.

# run them in reverse order, to be fair
restart:
with(LinearAlgebra):
with(Optimization):

N:=600:
U:=RandomMatrix(N,outputoptions=[datatype=float[8],shape=triangular[upper]]):
H:=Matrix(U.U^%T,shape=symmetric):

EV:=RandomVector(N,outputoptions=[datatype=float[8]]):

bl:=Vector(N,fill=0,datatype=float[8]):
bu:=Vector(N,fill=1,datatype=float[8]):

W:=Vector(N,symbol=w):
sol2 := CodeTools:-Usage( QPSolve(W^%T.H.W-W^%T.EV,seq(W[i]=bl[i]..bu[i],i=1..N)) ):
memory used=272.45MiB, alloc change=102.11MiB, cpu time=5.29s, real time=4.91s

sol1 := CodeTools:-Usage( QPSolve([-EV,2*H],[],[bl,bu]) ):
memory used=49.22MiB, alloc change=10.25MiB, cpu time=2.28s, real time=2.28s

Norm(Vector(sol1[2])-Vector(eval(W,sol2[2])));
                               0.

So Matrix form runs about twice as fast for this example, and seems to allocate only 70% more extra memory during the solving call. (Maple 15.01, 64bit Windows 7, intel i7).

acer

The short answer is that under `evalhf` coth is being evaluated for arguments outside of the range for which it computes accurately (without overflow, say). Internally, `plot` attempts with `evalhf` by default. You work around that problem when you wrap in `evalf` (or raise Digits above 15).

Examples of interest:

evalf(coth(355+I));

                                           -309  
               1.000000000 - 8.140551093 10     I

evalhf(coth(355+I));

                               1.

evalf(coth(356+I));

                                           -309  
               1.000000000 - 1.101703788 10     I

evalhf(coth(356+I)); # mildly interesting

             Float(undefined) + Float(undefined) I

evalf(1/(exp(354+I)));

                         -155                 -154  
           9.826304709 10     - 1.530356286 10     I

evalhf(1/(exp(354+I)));

                         -155                         -154  
   9.82630470844021329 10     - 1.53035628577376246 10     I

evalf(1/(exp(355+I)));

                         -155                 -155  
           3.614895484 10     - 5.629866151 10     I

evalhf(1/(exp(355+I))); # more interesting

                               0.

evalhf(1/(exp(355)*1/(exp(I))));

                         -155                         -155  
   3.61489548492129760 10     + 5.62986615203655732 10     I

acer

I'm not sure what you're intending for 'intercept', or if I've misunderstood you on that. I've assumed that you want to have 'intercept' be a parameter, where a single value for it is to be found. But if,instead, it's supposed to be yet another Vector of independent data then that can be accomodated below (if you supply it).

But Statistics:-Fit can handle your multiple independent data.

I guess you know that Vector GPA has only 48 entries, while the others have 50.

STARTING_SALARY := <29555, 27958, 27230, 31070, 27577, 30007, 28988, 32655, 29310,
 29145, 26142, 29460, 28390, 28574, 28760, 29665, 27481, 29186, 29302, 28917, 32839, 27277,
 26531, 32537, 28594, 28268, 29249, 28525, 26105, 30694, 33102, 30201, 26863, 31556, 29006,
 27054, 31677, 32430, 34023, 32786 , 34857, 34151, 35086, 33517, 36865, 33110, 32859,33181,
 34847, 34384>:

GPA := <2.74, 3.88, 2.70, 3.36, 2.57, 3.81, 2.70, 3.29, 2.97, 2.16, 2.84, 2.53, 2.79, 3.81,
 2.18, 2.25, 3.44, 3.04, 3.78, 2.03, 2.15, 3.89, 2.66, 3.07, 2.34, 3.89, 3.13, 3.40, 3.56, 2.74,
 2.04, 3.61, 2.62, 2.32, 2.35, 2.70, 3.89, 2.52, 3.35, 3.11, 3.72, 2.26, 2.23, 2.59, 2.19, 2.73,
 2.96, 2.66>:

METRICS := <0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0,
 0, 0, 0,0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1>:

SEX := <1, 0, 0, 1,1,0, 0, 0, 0, 1,0, 0, 1, 0, 0, 0, 0, 0, 1,1, 0, 0, 1,0, 0, 1, 0, 0, 1, 1,
 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1 ,1, 0, 0, 1>:

full := Statistics:-Fit( intercept + sex*B0 + B1*gpa + B2*metrics,
                 <GPA | METRICS[1..48] | SEX[1..48]>,
                 STARTING_SALARY[1..48],
                 [gpa, metrics, sex],
                 output=[leastsquaresfunction, parametervalues,
                         residualmeansquare, residuals] ):


form:=full[1]; # leastsquaresfunction

 59.004139262009176 sex - 694.6349921938285 gpa

    + 3978.586988002064 metrics

    + 31175.179243138406

full[2];

           [B0 = 59.004139262009176, 

             B1 = -694.6349921938285, 

             B2 = 3978.586988002064, 

             intercept = 31175.179243138406]

Vector(48, (i)->eval(form,
                     [sex=SEX[i],gpa=GPA[i],metrics=METRICS[i]])
                - STARTING_SALARY[i]): # this is just the residuals

LinearAlgebra:-Norm( full[4]^%T + % );

                                         -11
                   5.40012479177676141 10   

full[3]; # residualmeansquare

                                      6
                        3.838359311 10 

acer

This error message mostly occurs when one tries to assign to a formal parameter of a procedure, inside that procedure, when its value is not assignable.

Consider this example,

> restart;
> f := proc(x) x:=13; end proc:

> f([7]);
Error, (in f) illegal use of a formal parameter

What's happened is that Maple received a statement to assign 13 not to the name x but rather to the value of x. It was told to assign 13 to the list [7]. But that's not legal, as [7] is not assignable. Since x is a formal parameter of the procedure, that somewhat cryptic error message ensues. (...and it only makes sense to those whom it doesn't really help.)

There are a few ways to work around this. Usually, the attempt was a mistake in the first place, and a local variable should be used instead. If one really does want to assign to the name of the inbound argument (to replace its value, say) then this can be accomplished by passing in an uneval-quoted name, or by declaring that parameter as type uneval. Continuing the example above,

> z := [7]:

> f('z');

                               13

> z;

                               13

> #...or,

> f := proc(x::uneval) x:=13; end proc:

> z := [7]:

> f(z):

> z;

                               13

The above behaviour is sometimes described as the procedure having a side-effect on its argument (here, on name z).

There are a few other errors that are related to this case above,

> [7] := 13;
Error, illegal use of an object as a name

> 7 := 13;
Error, invalid left hand side of assignment

The first of those two has an entry in the Error Message Guide. All these errors should have such explanatory descriptions. They come up quite often. Yours just came up on Nov 21 on stackoverflow.

In my opinion, none of these error messages are extremely useful to the beginner. The error messages in the first of this last pair of example above and in your own example are unhelpful because they don't make it crystal clear that the "illegal use" is the attempted act of assignment. The error message of the second of the last examples above is unclear because it does not make crystal clear what is invalid about the left hand side of the assignment (ie. the fact that it is not assignable).

Banal as it might sound, perhaps better might be, "Error, assignment attempted to non-assignable object" for all three examples. That lets one know precisely what act was illegal, and what was invalid about the object.

Error message guide pages should show one the workarounds for the task at hand. If they have to explain the very wording and terms of the error message then the message is partly at fault.

acer

Are you saying that a function call has been assigned to the name `expr`?

> expr:=F(sin(a)-exp(b));
                          expr := F(sin(a) - exp(b))

> type(expr,function);   
                                     true

> op(0,expr);            
                                       F

> op(1,expr);            
                                sin(a) - exp(b)

Or, for an indexed function call,

> expr:=F[foo](sin(a)-exp(b));
                        expr := F[foo](sin(a) - exp(b))

> op(0,expr);                 
                                    F[foo]

> type(expr,function) and type(op(0,expr),'indexed');
                                     true

> op(0,op(0,expr));
                                       F

> op(op(0,expr));  
                                      foo

acer

If efficiency and speed matter here (ie. if the Matrix is large, or if this is critical code) then it's better to build up the results using successive multiplication than it is to generate all the powers separately.

If efficiency it's a concern for you, then you can sensibly ignore this.

Even using a smart "binary powering" technique for generating each higher power of an exact Matrix, it takes about 75-76 individual Matrix-Matrix multiplications in order to compute all powers from A^1 through to A^20. But of course it only takes 19 such individual multiplications in order to build them all up one at a time.

restart:

N:=10:
A:=LinearAlgebra:-RandomMatrix(N):

ba,st:=kernelopts(bytesalloc),time():
current:=LinearAlgebra:-IdentityMatrix(N):
for n from 1 to 20 do
   current:=A . current;
   #print('A'^n = current);
end do:
kernelopts(bytesalloc)-ba,time()-st;

                         1310480, 0.016

restart:

N:=10:
A:=LinearAlgebra:-RandomMatrix(N):

ba,st:=kernelopts(bytesalloc),time():
for n from 1 to 20 do
   current := A^n;
   #print('A'^n = current);
end do:
kernelopts(bytesalloc)-ba,time()-st;

                         2883056, 0.031

As N gets larger, the difference in costs may become more pronounced.

(sorry, folks, but I'm on a Horner's scheme kick right now...)

acer

One way to test equivalence of two expressions is to test whether their difference (subtraction) is equivalent (or even directly equal) to zero.

You could look at the help-pages evalb , is , simplify/details , testeq , and verify.

The command `Testzero` is underdocumented. By default `Testzero` uses `Normalizer`.

Some examples, which might be more clear after scanning those pages:

expr1:=sin(alpha+beta):
expr2:=sin(alpha)*cos(beta) + cos(alpha)*sin(beta):

evalb(simplify(expr1-expr2,'trig')=0);

                              true

Normalizer:=t->simplify(t,trig):

Testzero(expr1-expr2);

                              true

evalb(23423*2^5=32);

                             false

eqn:=sqrt(2)*sqrt(3)=6^(1/2);

                      (1/2)  (1/2)    (1/2)
                     2      3      = 6     

evalb(eqn);

                             false
is(eqn);

                              true

acer

See the documentation of the NAG Fortran Library function f08mef. This is like CLAPACK's dbdsqr function. It's documentation mentions the (QR based, or differential QD based) algorithms used.

You can see which NAG routine is used by LinearAlgebra, from userinfo messages.

> infolevel[LinearAlgebra]:=1:

> LinearAlgebra:-SingularValues(<>):
SingularValues: calling external function
SingularValues: NAG hw_f08kef
SingularValues: NAG hw_f08mef

Maple also ships with precompiled external libraries (shared objects clapack.dll or libclapack.so) containing the double-precision real and complex (dxxxx and zxxxx) routines from versions 3.0 of CLAPACK. If you want another routine from that library, you can call it yourself using the external-calling mechanism. (See here for an example of such user-defined wrapperless external-calling to CLAPACK from Maple.)

You can use these notes to help decide which CLAPACK function you'd like to use for singular value computation. Note that function dbdsdc for the divide and conquer algorithm is already avaliable in version 3.0. So you could call that directly from stock, shipped Maple, using define_external.

If you are interested in getting more state-of-the-art, then you could consider LAPACK 3.2.x or higher. See these release notes, and in particular section 2 for items 2.5 and 2.10 and the relevent references [5],[6]. Note that v.3.2.1 is avaliable aready as CLAPACK (which you can obtain as source or may be even precompiled). If you build such a more recent version of (C)LAPACK yourself (compiling, linking as shared dynamic library with appropriate exports, etc) then you would again be able to call its functions directly using the wrapperless forms of define_external.

acer