@Alec Mihailovs 

I do hope that it comes across that I am trying to be precise, and not rude. However, you didn't "just" show how to use evalhf with `eval`. You did imply at some point that it was evalhf which was causing the behaviour to be like that of double precision.

You explained, "evalhf produces the same result if the calculations were done in the hardware floats", apparently  as part of the rationale for wrapping it around `eval`. That implied that evalhf is what made it get the same result as double precision. I was correcting that erroneous implication (intended by you, or not).

In Maple 10, extracting a scalar from a float[8] rtable did not produce an HFloat scalar object, but it does in Maple 15. In Maple 10 such extraction produced a software float scalar with a similar kind of length as the doubles inside a float[8] rtable. And such software floats drop accuracy under arithmetic operations at lower working precision. And so the evalhf(eval(...)) approach fails in Maple 11, when `eval` became available to run under evalhf, as I showed by example.

As for Maple 15, the evalhf wrapping may not always be necessary under certain kinds of operation (eg. arithmetic) in order to get behaviour with higher accuracy (eg. like dbl prec.) that would be obtained by software floats at the same current Digits setting,

> # Maple 15. evalhf is not needed, to get behaviour like doubles on the NLPSolve result.
> restart: 

> evalf(subs(sol[2],1/3.));

                          0.3333333333

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                        [0., [x = 3.)]]

> eval(1/x,sol[2]);

                       0.3333333333333333

As for the bundled Maple Libraries, yes they have a mix, and no it is not necessarily wrong. I did say that evalf(...,n) was not wrong, used in the right way in the right hands. It's quite OK to do evalf(expr, n) where expr is assigned to be some expression and n is assigned to be some posint. It's clear that this will not hit the corner case I mentioned where only a preformed expression sequence is passed as single argument to evalf.

I earnestly hope that this has not come off to other readers as a personal attack by me. I've written on this site before, that I hold Alec's abilities -- Maple, coding, and math -- in high regard. And clearly there was a little confusion about the particularities of the poster's maple 10 version, so it seems good that various distinctions have been aired.

ps. "complete information"... is very amusing. Thanks for that.

@Alec Mihailovs One of my example's points is that your use of evalhf not relevant to how the result is obtained. The `eval` call within the `evalhf` call that you made in Maple 15 does not inherit the Digits setting (of evalhf).

As I mentioned, a result from `eval` similar to that using hardware doubles may be due to the NLPSolve result being an HFloat (and cascading of that through the arithmetic of the `eval`). But it's not because it is called from within evalhf.

For some subtle reasons (to do with a small amount of corner case examples, and so subtle that I forget them over and over again, and always have to comb email archives to re-find them...) the indexed form of `evalhf` is slightly superior to the two-argument calling sequence. I believe that's why the help page for evalf lists the former and not the latter in recent versions. It's understood, you have your own preference, and I for one will not criticize it.

[edited: ..had to dig, for this] Here's one such corner case. Some users, after getting in the habit of using the non-indexed form like evalf(expression, n) , get confused why these two calls below don't work the same way. It's not that evalf(..., n) is wrong, but more that people who aren't in the habit of using it won't run into this confusing example. (So, it may not affect an expert like yourself, but there's still a reasonable argument for suggesting one over the other to less experienced users.)

> evalf(1/3, 7);

                           0.3333333

> g:=1/3, 7;

                              1   
                              -, 7
                              3   

> evalf(g);

                        0.3333333333, 7.

I'm not a big fan of how votes change the order of posts here, either. Together with the mix and mix-ups of answers vs comments, it makes a mess of lots of threads.

@Alec Mihailovs One of my example's points is that your use of evalhf not relevant to how the result is obtained. The `eval` call within the `evalhf` call that you made in Maple 15 does not inherit the Digits setting (of evalhf).

As I mentioned, a result from `eval` similar to that using hardware doubles may be due to the NLPSolve result being an HFloat (and cascading of that through the arithmetic of the `eval`). But it's not because it is called from within evalhf.

For some subtle reasons (to do with a small amount of corner case examples, and so subtle that I forget them over and over again, and always have to comb email archives to re-find them...) the indexed form of `evalhf` is slightly superior to the two-argument calling sequence. I believe that's why the help page for evalf lists the former and not the latter in recent versions. It's understood, you have your own preference, and I for one will not criticize it.

[edited: ..had to dig, for this] Here's one such corner case. Some users, after getting in the habit of using the non-indexed form like evalf(expression, n) , get confused why these two calls below don't work the same way. It's not that evalf(..., n) is wrong, but more that people who aren't in the habit of using it won't run into this confusing example. (So, it may not affect an expert like yourself, but there's still a reasonable argument for suggesting one over the other to less experienced users.)

> evalf(1/3, 7);

                           0.3333333

> g:=1/3, 7;

                              1   
                              -, 7
                              3   

> evalf(g);

                        0.3333333333, 7.

I'm not a big fan of how votes change the order of posts here, either. Together with the mix and mix-ups of answers vs comments, it makes a mess of lots of threads.

@Alec Mihailovs I do not wish to give offense, but it looks to me as if the stated problem was with the precision of evaluation of the objective at the returned solution (as opposed to say the question of whether that solution were globally or locally optimal, etc).

In Maple 10, which the submitter mentioned using, NLPSolve does not return a result with HFloats. So not even double precision will be used under `eval`. Also, `eval` cannot be called from within `evalhf` in Maple 10.

In Maple 11, it still wouldn't be adequate, if greater working precision that the current Digits setting (eg, default of 10) were wanted:

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                               [0., [x = 3.]]

> evalf(subs(sol[2],1/x));

                                0.3333333333

> evalhf(eval(1/x,sol[2]));

                            0.333333333299999979

So, for the stated issue, in earlier Maple it's probably better to either raise Digits first. Or invoke evalhf with its indexed form, like evalf[n](...).

In Maple 15, double precision may get used, due to the presence of HFloats in the NLPSolve result.

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                    [0., [x = HFloat(3.0)]]

> evalf(subs(sol[2],1/x));

                   HFloat(0.3333333333333333)

> evalhf(eval(1/x,sol[2]));
                      0.333333333333333315

> lprint(sol);
[0., [x = HFloat(3.)]]

> kernelopts(version);

   Maple 15.00, X86 64 WINDOWS, Mar 21 2011, Build ID 582305

But if even greater working precision is required, then raising Digits above 15, or calling evalf[n](..) with n>15, may be required.

@Alec Mihailovs I do not wish to give offense, but it looks to me as if the stated problem was with the precision of evaluation of the objective at the returned solution (as opposed to say the question of whether that solution were globally or locally optimal, etc).

In Maple 10, which the submitter mentioned using, NLPSolve does not return a result with HFloats. So not even double precision will be used under `eval`. Also, `eval` cannot be called from within `evalhf` in Maple 10.

In Maple 11, it still wouldn't be adequate, if greater working precision that the current Digits setting (eg, default of 10) were wanted:

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                               [0., [x = 3.]]

> evalf(subs(sol[2],1/x));

                                0.3333333333

> evalhf(eval(1/x,sol[2]));

                            0.333333333299999979

So, for the stated issue, in earlier Maple it's probably better to either raise Digits first. Or invoke evalhf with its indexed form, like evalf[n](...).

In Maple 15, double precision may get used, due to the presence of HFloats in the NLPSolve result.

> A:=(x-3)^2:

> sol:=Optimization:-NLPSolve(A);

                    [0., [x = HFloat(3.0)]]

> evalf(subs(sol[2],1/x));

                   HFloat(0.3333333333333333)

> evalhf(eval(1/x,sol[2]));
                      0.333333333333333315

> lprint(sol);
[0., [x = HFloat(3.)]]

> kernelopts(version);

   Maple 15.00, X86 64 WINDOWS, Mar 21 2011, Build ID 582305

But if even greater working precision is required, then raising Digits above 15, or calling evalf[n](..) with n>15, may be required.

@Christopher2222 You do not have to have an existing (empty) worksheet saved, in order to open a new, empty one.

All you have to do is INTERFACE_WORKSHEET(':-display',':-text'=...) where the ... is a certain kind of long string of XML.

I don't think this avenue helps with Doug's question.

I think that this below shows that it is the kernel (not the GUI) which is preventing quit/done/stop from succeeding. It might be possible to bamboozle the kernel into thinking it was running in the commandline interface, but obvious attempts at that (like forcing a lie in the remember table of `IsWorksheetInterface`) don't seem to work here.

FromInert(_Inert_PROC(_Inert_PARAMSEQ(), _Inert_LOCALSEQ(),
                      _Inert_OPTIONSEQ(), _Inert_EXPSEQ(),
                      _Inert_STATSEQ(_Inert_STOP()),
                      _Inert_DESCRIPTIONSEQ(), _Inert_GLOBALSEQ(),
                      _Inert_LEXICALSEQ(), _Inert_EOP(_Inert_EXPSEQ())))();

Warning, done/quit/stop disabled.  Please use File->Close

What else could be tried? Hmm. Has anyone ever experimented with external-calling to Java on the libs bundled with the product? (Is there a public/private issue, with that? I'm also asking because I'd like to find the routine the GUI uses to encode embedded images -- it's not just base64.) I tried an OpenMaple help example once, and it acted as it it were CLI. Just wild thinking.

@Christopher2222 You do not have to have an existing (empty) worksheet saved, in order to open a new, empty one.

All you have to do is INTERFACE_WORKSHEET(':-display',':-text'=...) where the ... is a certain kind of long string of XML.

I don't think this avenue helps with Doug's question.

I think that this below shows that it is the kernel (not the GUI) which is preventing quit/done/stop from succeeding. It might be possible to bamboozle the kernel into thinking it was running in the commandline interface, but obvious attempts at that (like forcing a lie in the remember table of `IsWorksheetInterface`) don't seem to work here.

FromInert(_Inert_PROC(_Inert_PARAMSEQ(), _Inert_LOCALSEQ(),
                      _Inert_OPTIONSEQ(), _Inert_EXPSEQ(),
                      _Inert_STATSEQ(_Inert_STOP()),
                      _Inert_DESCRIPTIONSEQ(), _Inert_GLOBALSEQ(),
                      _Inert_LEXICALSEQ(), _Inert_EOP(_Inert_EXPSEQ())))();

Warning, done/quit/stop disabled.  Please use File->Close

What else could be tried? Hmm. Has anyone ever experimented with external-calling to Java on the libs bundled with the product? (Is there a public/private issue, with that? I'm also asking because I'd like to find the routine the GUI uses to encode embedded images -- it's not just base64.) I tried an OpenMaple help example once, and it acted as it it were CLI. Just wild thinking.

Thanks for adding the favorites. I am guessing that reading each others' favorites might become a popular pasttime.

I'm having trouble adding favorites to posts/answers which I've previously up-voted. (When I try, it just tells me that I've already voted, but doesn't add it as a favorite. It looks like interference with the up-vote icon right above the add-favorites icon.)

On my favorites page, any author's names seems to just be linked to www.maplesoft.com instead of to that user's profile page.

All the "Added" dates are not correct for anything submitted before "Primes v.2", ie. April 27, 2010.

Could the favourites please be sorted by the date of original submission of the respective posts?

The new Moderator badge is listed twice.

I'm sure that any issues can be fixed, so please don't take the above as criticism.

acer

@PatrickT Sure, it will work. But just like for other dsolve-numeric results with parameters, you'll have to set a parameter value before numerically computing.

Note the use of `forget`, since all of xt, yt, ft, and fft have been given option remember. This should be done each time you change any parameter value. If you don't clear their remember tables (by hand, or using `forget`) then their results will get remembered to correspond to the previous parameter value.

restart:

SOL := dsolve({ diff(x(t),t) = y(t)+x(t),
                diff(y(t),t) = y(t)-x(t),
                x(0) = p, y(0) = 1},
                'type' = numeric, 'output' = listprocedure,
                'parameters' = [p]):

xt := subsop(3=remember,eval( x(t), SOL )):
yt := subsop(3=remember,eval( y(t), SOL )):

f := (x,y) -> x^2+y^2:

ft := unapply( subs( X=xt(t), Y=yt(t), f(X,Y) ),
              t, 'proc_options' = ['remember','operator','arrow']):

temp := subs( unassigneddummy = f(xt(t),yt(t)), t -> evalf(unassigneddummy) ):
fft := subsop(3=op({op(3,eval(temp)),'remember'}), eval(temp) ):

ft(0.7);
Error, (in xt) parameters must be initialized before solution can be computed

fft(0.7);
Error, (in xt) parameters must be initialized before solution can be computed

SOL(parameters=[p=2.7]): #seems to work, but printed output looks odd

ft(0.7);

                              33.6176069516257

fft(0.7);

                              33.6176069516257

SOL(parameters=[p=3.4]):

ft(0.7); # rememebered!

                              33.6176069516257
fft(0.7); # remembered!

                              33.6176069516257


forget(ft):
forget(fft):
forget(xt):
forget(yt):

ft(0.7);

                              50.9333113513925
fft(0.7);

                              50.9333113513925

I wonder why doesn't searching for the words minimal and curvature find this interesting comment by DJKeenan.

acer

The posted code makes it easy to adjust the slopes at the knots. Thanks, Alec.

Let's recall that cubic spline interpolation, in a commonly used sense, minimizes curvature while producing a curve which is twice continuously differentiable. I guess people like this, if they hope to model smooth phenomena. There are lots of variants which can use piecewise cubics to produce something which may only be once continuously differentiable.

Christopher asked for the extra constraint that the slope be zero at certain points. Alec cleverly interpreted this to mean zero only at knots which are locally extreme (w.r.t the discrete adjacent knot values). This is a nice refinement on an instance of a Cardinal spline with parameter c=1, and also more sophisticated that what's called a Catmull–Rom (c=0) variant of the cubic Hermite spline.

Alec, I believe that your code does in fact use the average of the secants on either side, ie. the average of the slopes of the two line segments from knot i-1 to i and from knot i to i+1. (Your code gets this because the x-differences are all equal, I think... I tried to alter this to be true also for irregularly spaced knots, by using a scaled version of the so-called "finite difference" variant, see attached.)

But there's still a lot of leeway as to what tweaks one can make to the slope values at the knots. I'm attaching a worksheet in which the finite difference variant is applied to all internal knots which are not locally extreme. I keep Alec's twist to satisfy Chris, that slopes be zero at knots whose immediate neighbours are either both greater or lower in value. I scale each of the nonzero internal knots by parameter `c`. And then (unless I computed it wrongly) I find the `c` value which minimizes the total absolute curvature.

Here's a worksheet. The two animations should each produce the curve Alec showed above when c=1 (if I did it all right). Except that I changed the value at x=9 to be 7, instead of 6, so as to accentuate the hiccup that Christopher asked about in his comment above.

alec_modified_hermit.mw

Lots more variations are possible. One could use a set of c[i] and optimize all of them at once w.r.t the total absolute curvatue (each c[i] affects a pair of segments, just as the slopes themselves do). That might reduce hiccups and inflexions even more. One can think of even more interesting variants -- more than just those which have fancy sounding names.

I forgot to scale the slopes at the two end-knots by `c`, but changing the code is easy.

I didn't put any thought into how to compute the curvature more quickly (given that I know it's a piecewise, and I know the formula, etc). Corrections and comments are always welcome. The Optimization took about 10sec, and the plot a bit longer. There must be a fast way, but I didn't really try.

Alec's code layout is simple and easy to modify. I'm not sure whether the easier path to efficiency is to implement it using linear algebra, or to edit the procs to be Compilable.

acer

The spammed posts which have been cluttering up the Recent page make clear some other unfortunate aspects of the "new" Mapleprimes: a signigicant number of old posts have their formatting ruined, and are misclassified as Posts instead of Questions (and the replies as Comments instead of Answers).

The mis-format typically consists of each entire post or reply, including any pre-tag code, being crammed into one illegible paragraph.

I think that the main value of Mapleprimes is that it is a repository of great answers. But the broken search mechanism (indexing the "summary pages" as one problem has been called. Yikes!), the lack of good math typesetting, the difficulty in inserting 2D Math input except as attached full worksheets, the mis-format of old posts, the misclassification of old posts, and so on, take away from this value.

Mapleprimes could be more than just a free, one-stop (throw-away) Maple Support site, which is the epitome of a short-term benefit.

Fixing this site could transform it into an ongoing and easily accessed repository of Maple knowledge and solutions, and it could once again foster high quality technical discussion. Those are things which can bring additional long-term benefits such as increased web traffic, technical reputation, and good will.

Trying to think ahead: solving the problems don't (and really, really should not) need yet another re-implementation. Nobody wants another version 0.1 completely new site, twitter/facebook/cloud access, etc. Like for so many other things, the best course of action is to just amend what's already here with careful fixes to specific problems. Solid reliability using simple browser based access is better than version 0.1 implementation of new tech or even of old tech merged with new.

ps. I find that setting up Google Reader, with RSS feeds for the Question and Posts and Maplesoft Blog forums, is a convenient way to see all new postings without the spam getting too much in the way.

acer

@PatrickT 

You won't see the same order of speedup by giving Cru option remember as you might for xt and yt, unless you plan on calling Cru itself in several circumstances. But the methodology is pretty similar.

It doesn't arise here, but in general you might have to be careful about scoping, sometimes, and ensure that the names used in the unapply and funciton calls (eg. t) aren't assigned.

One interesting subtlety is that, if your original Cru is a procedure which does not apply evalf to its final result and if it contains exact numerics like sqrt(3), then you might want a resulting procedure which does in fact apply evalf. This doesn't affect plotting, as `plot` itself applies evalf or evalhf to final results. But you can see the distinction below between ct(0.7) and anotherct(0.7).

restart:

SOL:=dsolve({diff(x(t),t)=y(t)+x(t),diff(y(t),t)=y(t)-x(t),x(0)=1,y(0)=1},
             numeric,output=listprocedure):
xt := subsop(3=remember,eval( x(t), SOL )):
yt := subsop(3=remember,eval( y(t), SOL )):

# Let's suppose your original Cru is a procedure that takes two arguments (a,b)

Cru := (a,b) -> sqrt(3)/3*sqrt(1/(a+1/10))+sqrt(3)/90*sqrt(1/(b+1/10)^3):

ct := unapply(subs(A=xt(t),B=yt(t), Cru(A,B)),
              t, 'proc_options'=['remember','operator','arrow']):

lprint(eval(ct)); # note absence of evalf

ct(0.7);

temp := subs(unassigneddummy=Cru(xt(t),yt(t)), t->evalf(unassigneddummy)):

anotherct := subsop(3=op({op(3,eval(temp)),'remember'}), eval(temp) ):

lprint(eval(anotherct)); # note presence of evalf

anotherct(0.7);

plot(ct, 0.0 .. 0.7);
plot(anotherct, 0.0 .. 0.7);

@PatrickT You can place `option remember` on the procedures returned from dsolve/numeric. (Or you could wrap them in procedures that have option remember.)

Compare both times taken to produce the second plot, P2, below.

restart:
SOL:=dsolve({diff(x(t),t)=y(t)+x(t),diff(y(t),t)=y(t)-x(t),x(0)=1,y(0)=1},
             numeric,output=listprocedure):
xt := eval( x(t), SOL ):
yt := eval( y(t), SOL ):

f := (a,b)->sin(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P1:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=red ):
time()-st;

                                   0.080

f := (a,b)->cos(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P2:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=blue ):
time()-st;

                                   0.070

plots:-display(P1,P2);

restart:
SOL:=dsolve({diff(x(t),t)=y(t)+x(t),diff(y(t),t)=y(t)-x(t),x(0)=1,y(0)=1},
             numeric,output=listprocedure):
xt := subsop(3=remember,eval( x(t), SOL )):
yt := subsop(3=remember,eval( y(t), SOL )):

f := (a,b)->sin(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P1:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=red ):
time()-st;

                                   0.070

f := (a,b)->cos(a*b):
ft := T->f( xt(T), yt(T) ):

st:=time():
P2:=plots:-odeplot( SOL, [ ft(t), xt(t) ], color=blue ):
time()-st;

                                   0.010

plots:-display(P1,P2);

You can also force your plots to use specific domain points (ie. pointplot instead of odeplot) and thus enforce re-use of the same `t` values (and thus benefit from option remember) in subsequent plots of different `f` compositions.

Array output directly obtained from dsolve/numeric is a convenience. But you can also produce it yourself, after calling dsolve. So it might just be an inconvenience that dsolve & parameters doesn't allow ouput=Array, which you can work around.

This page, now spammed (June 08, 2011), deserves a screengrab.

acer