@Jimmy You can add 'parametervalues' to your 'output'=[...] and get the estimated parameter results, without having to solve for them using the final form.

@Jimmy You can add 'parametervalues' to your 'output'=[...] and get the estimated parameter results, without having to solve for them using the final form.

@herclau That doesn't prove much, since the original does not appear to be optimized for speed. Using a float[8] Array the procedure which populates the color values (inplace) could be Compiled or run under evalhf, after which it would be much faster.

Under evalhf, populating float[8] Array c inplace with a procedure, I see 3.4 times faster than with your approach.

Also, the original is not optimized for memory either. For N=25000 it can actually be done with only 3.3 MB allocated.

restart:
N:=25000:

xy:=LinearAlgebra:-RandomMatrix(N,2,generator=0.0..1.0,
                                outputoptions=[datatype=float[8]]):
str:=time[real]():
f:=proc(A::Array(datatype=float[8],order=C_order),
        xy::Matrix(datatype=float[8]), n)
 local i;
 for i from 1 to n*3 do
  A[i]:=`if`(irem(i,3)=1,xy[(i+2.0)/3.0,2],0);
 end do:
 NULL:
end proc:

p:=plots:-pointplot(xy,color=red,symbolsize=4):

c:=Array(1..N*3, datatype=float[8],order=C_order):
evalhf(f(c,xy,N)):
subsindets(p,specfunc(anything,COLOUR),z->'COLOUR'('RGB',c));

                           PLOT(..)  (actually displayed in my Maple 15)

time[real]()-str;

                             0.227

kernelopts(bytesalloc);

                            3276200

whereas

restart:
N:=25000:

xy:=LinearAlgebra:-RandomMatrix(N,2,generator=0.0..1.0,
                                outputoptions=[datatype=float[8]]):

str:=time[real]():

plots:-pointplot(xy,
                    color=[seq(RGB(xy[i, 2], 0., 0.), i = 1 .. N)],
                    symbolsize=4);


                           PLOT(..)  (actually displayed in my Maple 15)
time[real]()-str;

                              0.882

kernelopts(bytesalloc);
                            42721648

Using a single applicatiojn of `unapply` to create `g` can handle both Case 1 and 2 together.

> restart:

> mu := 0: sigma := 0.25: elist := [e1, e2, e3]:
> E := Statistics:-RandomVariable(Normal(mu, sigma)):
> jointPDF := product(Statistics:-PDF(E, elist[i]), i = 1..3):
> f := evalf(jointPDF):

> a := [e1+1.200000000 <= e2, e3-2.400000000 <= e1]:

> g := unapply(piecewise(a[1] and a[2], f, 0),[e1,e2,e3]): 

> evalf(Int( g, [-2 .. 2, -2 .. 2, -2 .. 2]), 3);
                                   0.000344

Using a single applicatiojn of `unapply` to create `g` can handle both Case 1 and 2 together.

> restart:

> mu := 0: sigma := 0.25: elist := [e1, e2, e3]:
> E := Statistics:-RandomVariable(Normal(mu, sigma)):
> jointPDF := product(Statistics:-PDF(E, elist[i]), i = 1..3):
> f := evalf(jointPDF):

> a := [e1+1.200000000 <= e2, e3-2.400000000 <= e1]:

> g := unapply(piecewise(a[1] and a[2], f, 0),[e1,e2,e3]): 

> evalf(Int( g, [-2 .. 2, -2 .. 2, -2 .. 2]), 3);
                                   0.000344

Using `unapply` can also allow this Case 1 to work as well.

> restart:                                                                     
> a := [e1+1.200000000 <= e2, e3-2.400000000 <= e1]:   
                        
> g := unapply(piecewise(a[1] and a[2],
>                 evalf(64.0*exp((-1)*8.0*e1^2)*exp((-1)    
>                       *8.0*e2^2)*exp((-1)*8.0*e3^2)/Pi^2),    
>                        0),
>              [e1,e2,e3]):

> evalf(Int( g, [-2 .. 2, -2 .. 2, -2 .. 2]), 3);  
                            
                                   0.000550

See now,

lprint(eval(g));
(e1, e2, e3) -> piecewise(e1+1.200000000 <= e2 and e3-2.400000000 <= e1,
6.484555750*exp(-8.0*e1^2)*exp(-8.0*e2^2)*exp(-8.0*e3^2),0)

Using `unapply` can also allow this Case 1 to work as well.

> restart:                                                                     
> a := [e1+1.200000000 <= e2, e3-2.400000000 <= e1]:   
                        
> g := unapply(piecewise(a[1] and a[2],
>                 evalf(64.0*exp((-1)*8.0*e1^2)*exp((-1)    
>                       *8.0*e2^2)*exp((-1)*8.0*e3^2)/Pi^2),    
>                        0),
>              [e1,e2,e3]):

> evalf(Int( g, [-2 .. 2, -2 .. 2, -2 .. 2]), 3);  
                            
                                   0.000550

See now,

lprint(eval(g));
(e1, e2, e3) -> piecewise(e1+1.200000000 <= e2 and e3-2.400000000 <= e1,
6.484555750*exp(-8.0*e1^2)*exp(-8.0*e2^2)*exp(-8.0*e3^2),0)

@PatrickT Are the 3rd and 4th images in your comment accidentally switched?

Why do you think that they should all come out the exactly the same?

@bernie Well, can't you just create a piecewise function, using whatever `Fit` returns only on the subrange for which you consider it valid? Elsewhere, you'd have the piecewise evaluate to something else.

@bernie Well, can't you just create a piecewise function, using whatever `Fit` returns only on the subrange for which you consider it valid? Elsewhere, you'd have the piecewise evaluate to something else.

Are you operating in 2D Math input mode? If so, is that an extra space between the solve and the opening round-bracket? If so, then remove that space, since the 2D Math parse will interpret solve (x-2=0, x) as solve*(x-2=0, x) where the extra space is denoting implicit multiplication.

Are you operating in 2D Math input mode? If so, is that an extra space between the solve and the opening round-bracket? If so, then remove that space, since the 2D Math parse will interpret solve (x-2=0, x) as solve*(x-2=0, x) where the extra space is denoting implicit multiplication.

Initially you claimed that this is an example of a "minor" change. And you show a one-line change in the `pdsolve` routine.

But then you claim that there must be additional, internal changes because the simple one-line edit is not enough to fix earlier releases. And you have no idea how extensive those additional changes might be. So this might not be a simple example after all.

Are you asking how to assemble a collection of FromInert/edit/ToInert changes, that would bring the Maple 12 code up to the functionality of Maple 16?

How would you disseminate such a patch, if it were based on a later version's code, without violating copyright?

How would you ensure that the patched Maple 12 would continue to properly run its own pdsolve, both on larger examples as well as on unit tests and unit failure checks? Why would it be prudent to load such a patch if it were not supported by a proper test run? A proper test run is not the example which was fixed. A proper test run is a battery that shows that the change had not broken usual desired behaviour against other examples.

I remembered that someone else asked this a while back, but it took me a while to find it.

In terms of the Help system (apart from the manuals) it's a little too difficult to uncover this stuff. The Help page for topic colondash isn't cross-referenced as much as it could be. The content and examples on that page are also mostly about package members' names. It'd be improved with another example explaining the subject of this thread, eg. for option keywords of procedures.