This sounds similar to a known issue with Maple 9.5, which was addressed in the point-release 9.5.2.

acer

Those advanced mathematical propositions are not within the scope of evalb's design or purpose.

They are not even within the current scope of the "mathematical" verifier `is`. (I don't see a continuous property or type, so I suppose that `is` couldn't make use of an assumption such as D(F)::continuous.) It would be lovely to learn that I'm wrong on that point.

acer

Those advanced mathematical propositions are not within the scope of evalb's design or purpose.

They are not even within the current scope of the "mathematical" verifier `is`. (I don't see a continuous property or type, so I suppose that `is` couldn't make use of an assumption such as D(F)::continuous.) It would be lovely to learn that I'm wrong on that point.

acer

The loop has n go from 1 to N.

Each time through the loop it computes y[n+1] for the current value of n.

The last (Nth) time through the loop, n=N.

So the last (Nth) time through, it computes y[N+1].

acer

The loop has n go from 1 to N.

Each time through the loop it computes y[n+1] for the current value of n.

The last (Nth) time through the loop, n=N.

So the last (Nth) time through, it computes y[N+1].

acer

I figured it out. Suppose that I'm editing a page as Filtered HTML, and even Preview it. Then I accidentally click and go to some other page. Then I hit the back-button to return to the editor & Preview. At that point, the Input format is switched to Plain text and I must remember to set it back to Filtered HTML. It's a minor bug.

Thanks very much for switching it for the above post, Scott.

acer

I am sure that I posted the above item with "input format" as "Filtered HTML", because I checked the links directly from the Previewed display. But now it appears (to me, at least) as plaintext.

acer

Yes, this is numerical DE solving.

Some links to help in understanding the relationship to numerical integration are here and here. Those pages have even more links to follow on the topic.

acer

Yes, this is numerical DE solving.

Some links to help in understanding the relationship to numerical integration are here and here. Those pages have even more links to follow on the topic.

acer

It appears as if he wants to use the Classic GUI. Since that's not available on OSX, he thought that perhaps he could use Maple for Linux instead.

As someone pointed out, Linux (ELF format) binaries don't run natively on an OSX machine (even if its architecture is Intel) so he might instead try running Linux as a virtual operating system.

I guess that there isn't any sort of mechanism to run Linux binaries directly under OSX, in say the same way that lxrun could run them on an Intel-Solaris system. Who knows, maybe lxrun will get ported to OSX...?

acer

There is a possible undesirable side-effect of using `assign` to zip together two Matrices.

Recall, the suggestion was zip(assign,A,B) or similar, where Matrix A had entries like a[1,1], etc. Suppose that a[1,1] were already assigned a symbolic (name) as its value. The zip would then assign to that underlying name! That's probably not what you might want to bring about.

For example, look what happens to name d below,

> A:=Matrix(2,2,symbol=a):
> d:='d':
> a[1,1]:=d:
> B:=Matrix([[3,0],[4,4]]):
> zip(assign,A,B):
> A, a[1,1], d;

                                [3    0]
                                [      ], 3, 3
                                [4    4]

But maybe we can construct a safer procedure to use for the zipping, which doesn't also assign to name d as a more extensive side-effect.

> restart:

> p := proc(a::evaln,b)
>   if type(eval(a,2),name) then
>     assign(eval(a,2),b);
>   end if:
>   NULL;
> end proc:

> A:=Matrix(2,2,symbol=a):
> d:='d':
> a[1,1]:=d:
> B:=Matrix([[3,0],[4,4]]):
> zip(p,A,B):
> A, a[1,1], d;
                                [3    0]
                                [      ], 3, d
                                [4    4]
 

Now what about your claim, that assignments into the named entries of C don't carry through to those of A. As long as the named entries are the same, then it does seem to work fine. For example,

> restart:

> p := proc(a::evaln,b)
>   if type(eval(a,2),name) then
>     assign(eval(a,2),b);
>   end if:
>   NULL;
> end proc:

> A:=Matrix(3,3,symbol=a):
> C:=Matrix(2,2,symbol=a):
> d:='d':
> a[1,1]:=d:
> B:=Matrix([[3,0],[4,4]]):
> zip(p,C,B):
> A, a[1,1], d;

                     [   3          0       a[1, 3]]
                     [                             ]
                     [   4          4       a[2, 3]], 3, d
                     [                             ]
                     [a[3, 1]    a[3, 2]    a[3, 3]]

All I can think of, is that maybe your C and A contained different symbolic 'a'. Maybe one of them contained the members of global table 'a' while the other had some local 'a', or maybe each had different local 'a' entries.

Now on to your example where the zipping of `assign` didn't work when target C had sparse storage.

> restart:

> C:=Matrix(2,2,[[a[1,1],0],[a[2,1],0]],storage=sparse);

                                   [a[1, 1]    0]
                              C := [            ]
                                   [a[2, 1]    0]
 
> B:=Matrix([[3,0],[4,4]]):

> zip(assign,C,B):
Error, (in assign) invalid arguments

> C;

                                [   3       0]
                                [            ]
                                [a[2, 1]    0]

The assignment to C failed after the (1,1) entry was assigned to. But it couldn't assign to 0, and aborted with an error. The examples above hint at why it might fail. It's unable to bring about the more extensive side-effect on 0 which is not a name. The procedure `p` helps here too.

> restart:

> C:=Matrix(2,2,[[a[1,1],0],[a[2,1],0]],storage=sparse);

                                   [a[1, 1]    0]
                              C := [            ]
                                   [a[2, 1]    0]

> B:=Matrix([[3,0],[4,4]]):

> p := proc(a::evaln,b)
>   if type(eval(a,2),name) then
>     assign(eval(a,2),b);
>   end if:
>   NULL;
> end proc:

> zip(p,C,B):
> C;

                                   [3    0]
                                   [      ]
                                   [4    0]

acer

There is a possible undesirable side-effect of using `assign` to zip together two Matrices.

Recall, the suggestion was zip(assign,A,B) or similar, where Matrix A had entries like a[1,1], etc. Suppose that a[1,1] were already assigned a symbolic (name) as its value. The zip would then assign to that underlying name! That's probably not what you might want to bring about.

For example, look what happens to name d below,

> A:=Matrix(2,2,symbol=a):
> d:='d':
> a[1,1]:=d:
> B:=Matrix([[3,0],[4,4]]):
> zip(assign,A,B):
> A, a[1,1], d;

                                [3    0]
                                [      ], 3, 3
                                [4    4]

But maybe we can construct a safer procedure to use for the zipping, which doesn't also assign to name d as a more extensive side-effect.

> restart:

> p := proc(a::evaln,b)
>   if type(eval(a,2),name) then
>     assign(eval(a,2),b);
>   end if:
>   NULL;
> end proc:

> A:=Matrix(2,2,symbol=a):
> d:='d':
> a[1,1]:=d:
> B:=Matrix([[3,0],[4,4]]):
> zip(p,A,B):
> A, a[1,1], d;
                                [3    0]
                                [      ], 3, d
                                [4    4]
 

Now what about your claim, that assignments into the named entries of C don't carry through to those of A. As long as the named entries are the same, then it does seem to work fine. For example,

> restart:

> p := proc(a::evaln,b)
>   if type(eval(a,2),name) then
>     assign(eval(a,2),b);
>   end if:
>   NULL;
> end proc:

> A:=Matrix(3,3,symbol=a):
> C:=Matrix(2,2,symbol=a):
> d:='d':
> a[1,1]:=d:
> B:=Matrix([[3,0],[4,4]]):
> zip(p,C,B):
> A, a[1,1], d;

                     [   3          0       a[1, 3]]
                     [                             ]
                     [   4          4       a[2, 3]], 3, d
                     [                             ]
                     [a[3, 1]    a[3, 2]    a[3, 3]]

All I can think of, is that maybe your C and A contained different symbolic 'a'. Maybe one of them contained the members of global table 'a' while the other had some local 'a', or maybe each had different local 'a' entries.

Now on to your example where the zipping of `assign` didn't work when target C had sparse storage.

> restart:

> C:=Matrix(2,2,[[a[1,1],0],[a[2,1],0]],storage=sparse);

                                   [a[1, 1]    0]
                              C := [            ]
                                   [a[2, 1]    0]
 
> B:=Matrix([[3,0],[4,4]]):

> zip(assign,C,B):
Error, (in assign) invalid arguments

> C;

                                [   3       0]
                                [            ]
                                [a[2, 1]    0]

The assignment to C failed after the (1,1) entry was assigned to. But it couldn't assign to 0, and aborted with an error. The examples above hint at why it might fail. It's unable to bring about the more extensive side-effect on 0 which is not a name. The procedure `p` helps here too.

> restart:

> C:=Matrix(2,2,[[a[1,1],0],[a[2,1],0]],storage=sparse);

                                   [a[1, 1]    0]
                              C := [            ]
                                   [a[2, 1]    0]

> B:=Matrix([[3,0],[4,4]]):

> p := proc(a::evaln,b)
>   if type(eval(a,2),name) then
>     assign(eval(a,2),b);
>   end if:
>   NULL;
> end proc:

> zip(p,C,B):
> C;

                                   [3    0]
                                   [      ]
                                   [4    0]

acer

I was unable to use the Editor (and add content to a comment) when using Netscape 7.2 (or mozilla, or Epiphany) on Linux.

But now, with firefox 2.0.0.11 for i686 Linux the Editor works fine.

Thank you for all the cool new features, Will.

acer

The first thing to note is that (sorry, Doug, to contradict your Second iteration hint) the Statistics:-Sample calls must be kept out of any inner loops. They are the performance killers. Now, in previous postings here and here on this topic the parameters of the distribution were static. In that case, it proved important to utilize a somewhat underdocumented feature of Statistics:-Sample. The feature is that, when supplied a distribution but not a sample size, the Sample routine returns an actual procedure. That procedure can itself be run, with sample size as argument, and that procedure avoids a lot of the overhead of Sample. The idea is to repeatedly call that procedure returned by Sample, and not repeatedly call Sample. For example, one might do it like so,


S := Statistics:-Sample(Statistics:-Distribution(...)); # outside loop


Then inside the loops call S(n) for a sample of size n. But in your example, the Poisson distribution itself is parametrized. Until tonight, I didn't realize that the method described above could also be used in such a case, Let's see what we get when we call Sample on that, without any sample size argument supplied to it, and with an indeterminate name as the parameter.

S := Statistics:-Sample(Statistics:-RandomVariable(Poisson(x)));
showstat(S);

Gosh. Look at what showstat displays. It shows that S is a procedure which
passes both sample size n and parameter x to an external function (some
precompiled routine). And look, S doesn't have x declared as a local, so inside S x could get a value picked up by Maple's lexical scoping ability. Normally, I'd just recommend using that S inside the loop, because it has so much less overhead than does Statistics:-Sample. It's a little tricky, because S uses x via lexical scoping. But it's possible to do this if x is not assigned a numeric value before S is created while x is assigned a numeric value before S(1) is invoked.

lpmrA := proc (lambda, phi, alpha, N)
option hfloat;
local x, y, Y, i, pr, kop;
  Y:=Array(0..N,datatype=integer[4]);
  # Must create pr before assigning to local x !
  pr :=
Statistics:-Sample(Statistics:-RandomVariable(Poisson(x)));
  x := 1; y := 1;
  for i to N do
    x := lambda+phi*x+alpha*(y-x);
    y := pr(1)[1]; # pr gets x by lexical scoping.
    Y[i] := trunc(y);
  end do;
  convert( Y, list);
end proc:

But we can cheat a little, and do even better. Look at what that procedure S does. It creates its own local procedure and then calls it immediately. So each time S gets called within our loops that local procedure will also get recreated. Gosh. We could just create that procedure ourselves, just once, outside our loops, and use it repeatedly. To do that it's necessary to call a few Statistics routines that aren't exports of the Statistics package. That can be achieved by temporarily unsetting the opaque nature of modules.

lpmrB := proc (lambda, phi, alpha, N)
option hfloat;
local x, y, Y, i, pr, kop;
  Y:=Array(0..N,datatype=integer[4]);
  x := 1; y := 1;
  kop := kernelopts(opaquemodules);
  kernelopts(opaquemodules=false);
  Statistics:-ExternalSupport:-Initialize();
  pr :=
Statistics:-ExternalSupport:-DefineExternal("MaplePoissonRandomSample");
 kernelopts(opaquemodules=kop);
  for i to N do
    x := lambda+phi*x+alpha*(y-x);
    y := pr(1,x)[1];
    Y[i] := trunc(y);
  end do;
  convert( Y, list);
  end proc:

For me, lpmrA took about 0.16 sec while lpmrB took about half that time, to generate 3000 elements.

acer

The first thing to note is that (sorry, Doug, to contradict your Second iteration hint) the Statistics:-Sample calls must be kept out of any inner loops. They are the performance killers. Now, in previous postings here and here on this topic the parameters of the distribution were static. In that case, it proved important to utilize a somewhat underdocumented feature of Statistics:-Sample. The feature is that, when supplied a distribution but not a sample size, the Sample routine returns an actual procedure. That procedure can itself be run, with sample size as argument, and that procedure avoids a lot of the overhead of Sample. The idea is to repeatedly call that procedure returned by Sample, and not repeatedly call Sample. For example, one might do it like so,


S := Statistics:-Sample(Statistics:-Distribution(...)); # outside loop


Then inside the loops call S(n) for a sample of size n. But in your example, the Poisson distribution itself is parametrized. Until tonight, I didn't realize that the method described above could also be used in such a case, Let's see what we get when we call Sample on that, without any sample size argument supplied to it, and with an indeterminate name as the parameter.

S := Statistics:-Sample(Statistics:-RandomVariable(Poisson(x)));
showstat(S);

Gosh. Look at what showstat displays. It shows that S is a procedure which
passes both sample size n and parameter x to an external function (some
precompiled routine). And look, S doesn't have x declared as a local, so inside S x could get a value picked up by Maple's lexical scoping ability. Normally, I'd just recommend using that S inside the loop, because it has so much less overhead than does Statistics:-Sample. It's a little tricky, because S uses x via lexical scoping. But it's possible to do this if x is not assigned a numeric value before S is created while x is assigned a numeric value before S(1) is invoked.

lpmrA := proc (lambda, phi, alpha, N)
option hfloat;
local x, y, Y, i, pr, kop;
  Y:=Array(0..N,datatype=integer[4]);
  # Must create pr before assigning to local x !
  pr :=
Statistics:-Sample(Statistics:-RandomVariable(Poisson(x)));
  x := 1; y := 1;
  for i to N do
    x := lambda+phi*x+alpha*(y-x);
    y := pr(1)[1]; # pr gets x by lexical scoping.
    Y[i] := trunc(y);
  end do;
  convert( Y, list);
end proc:

But we can cheat a little, and do even better. Look at what that procedure S does. It creates its own local procedure and then calls it immediately. So each time S gets called within our loops that local procedure will also get recreated. Gosh. We could just create that procedure ourselves, just once, outside our loops, and use it repeatedly. To do that it's necessary to call a few Statistics routines that aren't exports of the Statistics package. That can be achieved by temporarily unsetting the opaque nature of modules.

lpmrB := proc (lambda, phi, alpha, N)
option hfloat;
local x, y, Y, i, pr, kop;
  Y:=Array(0..N,datatype=integer[4]);
  x := 1; y := 1;
  kop := kernelopts(opaquemodules);
  kernelopts(opaquemodules=false);
  Statistics:-ExternalSupport:-Initialize();
  pr :=
Statistics:-ExternalSupport:-DefineExternal("MaplePoissonRandomSample");
 kernelopts(opaquemodules=kop);
  for i to N do
    x := lambda+phi*x+alpha*(y-x);
    y := pr(1,x)[1];
    Y[i] := trunc(y);
  end do;
  convert( Y, list);
  end proc:

For me, lpmrA took about 0.16 sec while lpmrB took about half that time, to generate 3000 elements.

acer