The objective may not actually be the sum of the distances, but the command as written does provide a result which attains the minimum of that sum, no? That was my intended point, by saying that the code "minimized" that sum.

I would also suggest that its objective is less complicated than that alternative add() call, by virtue of not containing a call to sqrt.

acer

The objective may not actually be the sum of the distances, but the command as written does provide a result which attains the minimum of that sum, no? That was my intended point, by saying that the code "minimized" that sum.

I would also suggest that its objective is less complicated than that alternative add() call, by virtue of not containing a call to sqrt.

acer

The commands that I wrote had to all be done in exacty the order that I gave them. I suspect that you repeated the intersection() command.

Many of the geom3d commands work by assigning to a name (parameter) as a "side effect". If you do try to repeat the same command then that name gets evaluated to its previously assigned value when passed again as an argument to the procedure. Ie, the routine sees the previously assigned value coming in as the argument, and not an actual name as expected. Maple complains about this with the error message that you mentioned. You can work around this by simply quoting the relevant arguments. I give an example below.

> with(geom3d):
> point(a,-3,-10,2),point(b,5,14,8):
> line(L,[a,b]):
> _EnvXName := x: _EnvYName := y: _EnvZName := z:
> sphere(s,x^2+y^2+z^2 = 36,[x,y,z]):

> intersection(P,s,L): # this assigns to P
areinterls:   "two points of intersection"

> intersection(P,s,L): # it sees the value of P, not the name P
areinterls:   "two points of intersection"
Error, (in geom3d:-intersection) illegal use of a formal parameter

> intersection('P',s,L): # so pass in the name P
areinterls:   "two points of intersection"

You also can quote the relevant names the first (and every) time the command gets called, so that the code for it always looks similar.

The `areinterls: "two points of intersection"' bit is just a printed informative message, not a returned result. It's letting you know that it found two distinct points of intersection, that's all.

In order to understand how I showed it being done with Maple, it would be best to understand mathematically and geometrically the explanation I gave first about how to do it by hand.

acer

The commands that I wrote had to all be done in exacty the order that I gave them. I suspect that you repeated the intersection() command.

Many of the geom3d commands work by assigning to a name (parameter) as a "side effect". If you do try to repeat the same command then that name gets evaluated to its previously assigned value when passed again as an argument to the procedure. Ie, the routine sees the previously assigned value coming in as the argument, and not an actual name as expected. Maple complains about this with the error message that you mentioned. You can work around this by simply quoting the relevant arguments. I give an example below.

> with(geom3d):
> point(a,-3,-10,2),point(b,5,14,8):
> line(L,[a,b]):
> _EnvXName := x: _EnvYName := y: _EnvZName := z:
> sphere(s,x^2+y^2+z^2 = 36,[x,y,z]):

> intersection(P,s,L): # this assigns to P
areinterls:   "two points of intersection"

> intersection(P,s,L): # it sees the value of P, not the name P
areinterls:   "two points of intersection"
Error, (in geom3d:-intersection) illegal use of a formal parameter

> intersection('P',s,L): # so pass in the name P
areinterls:   "two points of intersection"

You also can quote the relevant names the first (and every) time the command gets called, so that the code for it always looks similar.

The `areinterls: "two points of intersection"' bit is just a printed informative message, not a returned result. It's letting you know that it found two distinct points of intersection, that's all.

In order to understand how I showed it being done with Maple, it would be best to understand mathematically and geometrically the explanation I gave first about how to do it by hand.

acer

You might find this syntax to be straightforward,

[seq(x[1..3], x in L)];

acer

You might find this syntax to be straightforward,

[seq(x[1..3], x in L)];

acer

I've edited out this comment, which I realized wasn't right.

acer

I've edited out this comment, which I realized wasn't right.

acer

The plotter is being "smart" and ignoring the (relatively) small imaginary components generated when evaluating f_inverse(y).

Consider the plot  of the imaginary components alone.

Fiy := unapply(f_inverse(y),y):
plot(Im@Fiy,0..15,y=0..2e-10);

evalf[10](Fiy(4.0));

Have a look at what happens if you instead create it as,

> ef_inverse:= y -> Invert(2,13-8*x+4*x^2-x^3,3,
>      -97848849/2500000+40116283/625000*x
>      -290813981/10000000*x^2+40116283/10000000*x^3,
>      9,2003409407/250000000-173624697/50000000*x
>      +3065739839/10000000000*x^2-658003077/50000000000*x^3):

> eFiy:=unapply(simplify(ef_inverse(y)),y):

Then you could look at,

plot(evalf@Im@eFiy,0..5,y=0..2e-10);

Digits:=100:
plot(evalf@Im@eFiy,0..5,y=0..2e-10);

evalf[100](eFiy(4.0));

Alternatively, you could raise Digits before creating f_inverse(y), and then notice that the imaginary artefacts become smaller as a consequence.

acer

First the original, with 10000 repetitions,

> restart: kernelopts(printbytes=false):
> with(Statistics):
> st:=time():
> for i from 1 to 10000 do
> B:=2*(Sample(Bernoulli(0.5),2)-<0.5,0.5>):
> end do:
> time()-st;
                                    13.709

And now a modification, to pull some Statistics overhead outside the loop.

> restart: kernelopts(printbytes=false):
> with(Statistics):
> with(LinearAlgebra):
> S := Sample(Bernoulli(0.5)):
> V := Vector[row](2,[0.5,0.5],'datatype'='float[8]'):
> st:=time():
> for i from 1 to 10000 do
> B:=S(2);
> VectorScalarMultiply(VectorAdd(B,V,1,-1,'inplace'=true),2.0,'inplace'=true):
> end do:
> time()-st;
                                     3.684

What if we create all 10000 sample pairs at once?

> restart: kernelopts(printbytes=false):
> with(Statistics):
> with(LinearAlgebra):
> S := Sample(Bernoulli(0.5)):
> st:=time():
> B:=ArrayTools:-Alias(S(20000),[10000,2]):
> V:=Matrix(10000,2,fill=0.5,'datatype'='float[8]'):
> MatrixScalarMultiply(MatrixAdd(B,V,1,-1,'inplace'=true),2.0,'inplace'=true):
> time()-st;
                                     0.025

And now to produce that plot,

> restart;
> F:=proc(n)
> option hfloat;
> local x,S,B,V,k,nn;
> uses Statistics, LinearAlgebra;
> nn:=trunc(n);
> x:=Array(0..nn,'datatype'='float[8]');
> x[1]:= 1:
> S := Sample(Bernoulli(0.5)):
> B:=ArrayTools:-Alias(S(2*nn),[nn,2]);
> V:=Matrix(nn,2,fill=0.5,'datatype'='float[8]'):
> MatrixScalarMultiply(MatrixAdd(B,V,1,-1,'inplace'=true),2.0,'inplace'=true):
> for k from 2 to nn do
> x[k]:=B[k,1]*x[k-1]+B[k,2]*x[k-2];
> end do;
> end proc:

> S:= n-> abs(F(n))^(1/n):

> time( assign('p',plot(S,1..10)) );
                                     0.936

That same plot took 5.2 sec to create using the original code.

There are some other tricks, but that covers some of the most "re-usable" ones.

acer

First the original, with 10000 repetitions,

> restart: kernelopts(printbytes=false):
> with(Statistics):
> st:=time():
> for i from 1 to 10000 do
> B:=2*(Sample(Bernoulli(0.5),2)-<0.5,0.5>):
> end do:
> time()-st;
                                    13.709

And now a modification, to pull some Statistics overhead outside the loop.

> restart: kernelopts(printbytes=false):
> with(Statistics):
> with(LinearAlgebra):
> S := Sample(Bernoulli(0.5)):
> V := Vector[row](2,[0.5,0.5],'datatype'='float[8]'):
> st:=time():
> for i from 1 to 10000 do
> B:=S(2);
> VectorScalarMultiply(VectorAdd(B,V,1,-1,'inplace'=true),2.0,'inplace'=true):
> end do:
> time()-st;
                                     3.684

What if we create all 10000 sample pairs at once?

> restart: kernelopts(printbytes=false):
> with(Statistics):
> with(LinearAlgebra):
> S := Sample(Bernoulli(0.5)):
> st:=time():
> B:=ArrayTools:-Alias(S(20000),[10000,2]):
> V:=Matrix(10000,2,fill=0.5,'datatype'='float[8]'):
> MatrixScalarMultiply(MatrixAdd(B,V,1,-1,'inplace'=true),2.0,'inplace'=true):
> time()-st;
                                     0.025

And now to produce that plot,

> restart;
> F:=proc(n)
> option hfloat;
> local x,S,B,V,k,nn;
> uses Statistics, LinearAlgebra;
> nn:=trunc(n);
> x:=Array(0..nn,'datatype'='float[8]');
> x[1]:= 1:
> S := Sample(Bernoulli(0.5)):
> B:=ArrayTools:-Alias(S(2*nn),[nn,2]);
> V:=Matrix(nn,2,fill=0.5,'datatype'='float[8]'):
> MatrixScalarMultiply(MatrixAdd(B,V,1,-1,'inplace'=true),2.0,'inplace'=true):
> for k from 2 to nn do
> x[k]:=B[k,1]*x[k-1]+B[k,2]*x[k-2];
> end do;
> end proc:

> S:= n-> abs(F(n))^(1/n):

> time( assign('p',plot(S,1..10)) );
                                     0.936

That same plot took 5.2 sec to create using the original code.

There are some other tricks, but that covers some of the most "re-usable" ones.

acer

It looks like the same problem recurred, for your posted code. The same bit is missing.

You might consider toggling the Source button (in the menu right above the editor box, when you post), and then pasting the code in between <pre></pre> tags. And then right after you paste, and before toggling Source back or previewing or anything, changeinstances of < to &lt; all by hand.

Or just expand the "Input format" thingy, right below the input box, and toggle it as "Plain Text". And then enter it all.

acer

It looks like the same problem recurred, for your posted code. The same bit is missing.

You might consider toggling the Source button (in the menu right above the editor box, when you post), and then pasting the code in between <pre></pre> tags. And then right after you paste, and before toggling Source back or previewing or anything, changeinstances of < to &lt; all by hand.

Or just expand the "Input format" thingy, right below the input box, and toggle it as "Plain Text". And then enter it all.

acer

I'm glad that you've been able to make progress.

Now that you have found your own way, and for your amusement,

> restart:
> a:=proc(n) option remember;2+ln(procname(n-1));end proc:
> a(0):=1.0:
> for i while abs(a(i)-a(i-1))>=1.0e-5 do end do:
> i, a(i);
                                13, 3.146191522

acer

I'm glad that you've been able to make progress.

Now that you have found your own way, and for your amusement,

> restart:
> a:=proc(n) option remember;2+ln(procname(n-1));end proc:
> a(0):=1.0:
> for i while abs(a(i)-a(i-1))>=1.0e-5 do end do:
> i, a(i);
                                13, 3.146191522

acer