pagan

5137 Reputation

23 Badges

16 years, 229 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are replies submitted by pagan

Also, if the escaped local isn't due to a mistake (ie. was intentional..) then you can still test against x^2 by converting to global.

eg, if the expression with x*x is assigned to name `ee`, then you can call

convert(ee,`global`);

Since the way I proposed (essentially, like Keep1 here) and the way Robert proprosed (essentially, like Keep2 here) are pretty similar in resource use and simplicity, perhaps it's the predicate which could be the focus.

> restart:
> A := Matrix(1..20000,1..3,rand(-100..100)):

> Keep1 := M ->
            M[[seq(`if`(M[i,1] >= 0,i,NULL), i=1..op([1,1],M))],..]:

> CodeTools:-Usage( Keep1(A) ):
memory used=0.99MiB, alloc change=1.12MiB, cpu time=15.00ms, real time=15.00ms

> restart:
> A := Matrix(1..20000,1..3,rand(-100..100)):

> Keep2 := M ->
            M[remove(t -> (0 > M[t,1]), [$1..op([1,1],M)]),..]:

> CodeTools:-Usage( Keep2(A) ):
memory used=2.22MiB, alloc change=2.25MiB, cpu time=31.00ms, real time=32.00ms

Using `remove` costs a little more, above. But it looks easier to adjust Keep2 so that it would accept a second argument which could be a predicate, passed as an operator which is handy for `remove`.

nb. I'm running this now on 64bit Maple, which explains why Keep1 is using twice the memory as my earlier MatrixSelectRowsA which I ran on 32bit Maple.

Test results (together). Time them separately, else the performance (esp. w.r.t. memory) is often affected by which is done first, etc.

restart:
M := Matrix(1..200,1..3,rand(-100..100)):
MatrixSelectRows1 := proc(M::Matrix)
  local i ;
  Matrix(=0,M[i,1..-1],NULL),i=1..op([1,1],M))>);
end proc:
M21 := MatrixSelectRows1(M):
MatrixSelectRows2 := proc(M::Matrix)
  local i, U ;
  U := NULL :
  for i to LinearAlgebra:-RowDimension(M) do
    if M[i,1] < 0
      then U := U,i ;
    end if ;
  end do ;
  LinearAlgebra:-DeleteRow(M,[U]);
end proc:
M22 := MatrixSelectRows2(M):
LinearAlgebra:-Norm(M22-M21);

                               0

MatrixSelectRowsA := proc(M::Matrix)
  local i, U, Ucompl, m ;
  m := op([1,1],M);
  U := seq(`if`(M[i,1]>=0,i,NULL),i=1..m);
  M[[U],1..-1];
end proc:
M2A := MatrixSelectRowsA(M):
LinearAlgebra:-Norm(M2A-M21);

                               0

restart:
M := Matrix(1..20000,1..3,rand(-100..100)):
MatrixSelectRows1 := proc(M::Matrix)
  local i ;
  Matrix(=0,M[i,1..-1],NULL),i=1..op([1,1],M))>);
end proc:
st,ba,bu := time(),kernelopts(bytesalloc),kernelopts(bytesused):
M21 := MatrixSelectRows1(M):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

                    0.063, 1965720, 2174084

restart:
MatrixSelectRows2 := proc(M::Matrix)
  local i, U ;
  U := NULL :
  for i to LinearAlgebra:-RowDimension(M) do
    if M[i,1] < 0
      then U := U,i ;
    end if ;
  end do ;
  LinearAlgebra:-DeleteRow(M,[U]);
end proc:
M := Matrix(1..20000,1..3,rand(-100..100)):
st,ba,bu := time(),kernelopts(bytesalloc),kernelopts(bytesused):
M22 := MatrixSelectRows2(M):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

                   0.515, 70569348, 200748624

restart:
MatrixSelectRowsA := proc(M::Matrix)
  local i, U, Ucompl, m ;
  m := op([1,1],M);
  U := seq(`if`(M[i,1]>=0,i,NULL),i=1..m);
  M[[U],1..-1];
end proc:
M := Matrix(1..20000,1..3,rand(-100..100)):
st,ba,bu := time(),kernelopts(bytesalloc),kernelopts(bytesused):
M2A := MatrixSelectRowsA(M):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

                     0.031, 655240, 602196
eq1 := (a+1)*x - b*y + c*z = j:
eq2 := 2*x  + (sqrt(d)+e)*y - f*z = k:
eq3 := (g-h*2+1)*x - i*y - 4*z = l:

A,B := LinearAlgebra:-GenerateMatrix([eq1,eq2,eq3],[x,y,z]);

                         [   a + 1         -b      c ]  [j]
                         [                           ]  [ ]
                         [              (1/2)        ]  [k]
                 A, B := [     2       d      + e  -f], [ ]
                         [                           ]  [l]
                         [g - 2 h + 1      -i      -4] 
    
A.Vector([x,y,z])=B;

                     [   (a + 1) x - b y + c z   ]   [j]
                     [                           ]   [ ]
                     [      / (1/2)    \         ]   [k]
                     [2 x + \d      + e/ y - f z ] = [ ]
                     [                           ]   [l]
                     [(g - 2 h + 1) x - i y - 4 z] 

A.Vector([x,y,z])=~B;

                      [   (a + 1) x - b y + c z = j   ]
                      [                               ]
                      [      / (1/2)    \             ]
                      [2 x + \d      + e/ y - f z = k ]
                      [                               ]
                      [(g - 2 h + 1) x - i y - 4 z = l]
     
LinearAlgebra:-GenerateEquations(A,[x,y,z],B):

map(print,%):

                           x a + x - b y + c z = j
                                (1/2)                
                       2 x + y d      + y e - f z = k
                       x g - 2 x h + x - i y - 4 z = l
eq1 := (a+1)*x - b*y + c*z = j:
eq2 := 2*x  + (sqrt(d)+e)*y - f*z = k:
eq3 := (g-h*2+1)*x - i*y - 4*z = l:

A,B := LinearAlgebra:-GenerateMatrix([eq1,eq2,eq3],[x,y,z]);

                         [   a + 1         -b      c ]  [j]
                         [                           ]  [ ]
                         [              (1/2)        ]  [k]
                 A, B := [     2       d      + e  -f], [ ]
                         [                           ]  [l]
                         [g - 2 h + 1      -i      -4] 
    
A.Vector([x,y,z])=B;

                     [   (a + 1) x - b y + c z   ]   [j]
                     [                           ]   [ ]
                     [      / (1/2)    \         ]   [k]
                     [2 x + \d      + e/ y - f z ] = [ ]
                     [                           ]   [l]
                     [(g - 2 h + 1) x - i y - 4 z] 

A.Vector([x,y,z])=~B;

                      [   (a + 1) x - b y + c z = j   ]
                      [                               ]
                      [      / (1/2)    \             ]
                      [2 x + \d      + e/ y - f z = k ]
                      [                               ]
                      [(g - 2 h + 1) x - i y - 4 z = l]
     
LinearAlgebra:-GenerateEquations(A,[x,y,z],B):

map(print,%):

                           x a + x - b y + c z = j
                                (1/2)                
                       2 x + y d      + y e - f z = k
                       x g - 2 x h + x - i y - 4 z = l

It's a small data set of only three points, but I wondered whether he was hoping for a best fit of the data to the curve a*(x-c)^b=y

It's a small data set of only three points, but I wondered whether he was hoping for a best fit of the data to the curve a*(x-c)^b=y

That's useful.

Sometimes it might be desirable to restrict the expanding.

dividend2:=proc(ex,f:=numer(ex))
local n,d;
    n:=frontend(expand,[numer(ex)/f]):
    d:=frontend(expand,[denom(ex)/f]);
    n/d;
end proc:

A simple example to compare might be (R*D^2)/(R*D^2+s*cos(a+b+c)+s^2*R*exp(z-t)*C)

That's useful.

Sometimes it might be desirable to restrict the expanding.

dividend2:=proc(ex,f:=numer(ex))
local n,d;
    n:=frontend(expand,[numer(ex)/f]):
    d:=frontend(expand,[denom(ex)/f]);
    n/d;
end proc:

A simple example to compare might be (R*D^2)/(R*D^2+s*cos(a+b+c)+s^2*R*exp(z-t)*C)

@PatrickT That Mandelbrot Mania post has code which is too slow.

That "Mania" worksheet uses `option hfloat` to try and gain speed. Yet that worksheet even contains this sentence, "Please note that this procedure can take some time to execute; you may want to start with a low value of pts (perhaps 300) before creating larger pictures."

Using the faster `evalhf` floating-point interpreter on an procedure which writes to the float[8] Array "in-place" outperforms `option hfloat` significantly here. That generates no garbage. Option hfloat causes HFloat objects to be created, and those (each) are DAGS and as such collectible garbage, which is one reason it's slow. (Other reasons are that all other flow control is still slow...)

Much faster still is an in-place procedure used with the Maple Compiler. See here for such an approach. I forget the exact details but if I recall this is something like up to a 1000 times faster than that "Mania" worksheet's code on the main example.

It varies with the application, but sometimes evalhf can be about 10-15 times faster than regular Maple (including when using `option hfloat`, as here), and the Compiler can sometimes bring another 10-20 times speed up again. And both can be used effectiely to reduce the expense of memory use.

 

@PatrickT That Mandelbrot Mania post has code which is too slow.

That "Mania" worksheet uses `option hfloat` to try and gain speed. Yet that worksheet even contains this sentence, "Please note that this procedure can take some time to execute; you may want to start with a low value of pts (perhaps 300) before creating larger pictures."

Using the faster `evalhf` floating-point interpreter on an procedure which writes to the float[8] Array "in-place" outperforms `option hfloat` significantly here. That generates no garbage. Option hfloat causes HFloat objects to be created, and those (each) are DAGS and as such collectible garbage, which is one reason it's slow. (Other reasons are that all other flow control is still slow...)

Much faster still is an in-place procedure used with the Maple Compiler. See here for such an approach. I forget the exact details but if I recall this is something like up to a 1000 times faster than that "Mania" worksheet's code on the main example.

It varies with the application, but sometimes evalhf can be about 10-15 times faster than regular Maple (including when using `option hfloat`, as here), and the Compiler can sometimes bring another 10-20 times speed up again. And both can be used effectiely to reduce the expense of memory use.

 

@PatrickT One of the things to notice about this answer involving % is that a procedure U is only created once. This is one of the things that makes it a good answer.

Any other solution which involves creating new instances of a procedure U every time you want to change which parameter values are used, is bad in some respects.

Notice that for this answer the procedure is formed once and is assigned directly to U, which is also true in at least one other answer.

Also, the piecewise solution doesn't look shorter or simpler here either:

Functions2 := [ U = 'proc(c) global s; piecewise(s=1,log(c),(c^(1-s)-1)/(1-s)) end proc' ]:

Parameters2 := [ s = 1 ];
                            [s = 1]

eval(eval(U(c), Functions2), Parameters2);
                             ln(c)

Parameters2 := [ s = 2 ];
                            [s = 2]

eval(eval(U(c), Functions2), Parameters2);
                              1    
                            - - + 1
                              c    

# Compare with..

U := proc(C, S) if S=1 then log(C); else (C^(1-S)-1)/(1-S); end if; end proc:

Parameters1 := [ s = 1 ];
                            [s = 1]

U(c, subs(Parameters1,s));
                             ln(c)

Parameters1 := [ s = 2 ];
                            [s = 2]

U(c, subs(Parameters1,s));
                              1    
                            - - + 1
                              c    

@PatrickT One of the things to notice about this answer involving % is that a procedure U is only created once. This is one of the things that makes it a good answer.

Any other solution which involves creating new instances of a procedure U every time you want to change which parameter values are used, is bad in some respects.

Notice that for this answer the procedure is formed once and is assigned directly to U, which is also true in at least one other answer.

Also, the piecewise solution doesn't look shorter or simpler here either:

Functions2 := [ U = 'proc(c) global s; piecewise(s=1,log(c),(c^(1-s)-1)/(1-s)) end proc' ]:

Parameters2 := [ s = 1 ];
                            [s = 1]

eval(eval(U(c), Functions2), Parameters2);
                             ln(c)

Parameters2 := [ s = 2 ];
                            [s = 2]

eval(eval(U(c), Functions2), Parameters2);
                              1    
                            - - + 1
                              c    

# Compare with..

U := proc(C, S) if S=1 then log(C); else (C^(1-S)-1)/(1-S); end if; end proc:

Parameters1 := [ s = 1 ];
                            [s = 1]

U(c, subs(Parameters1,s));
                             ln(c)

Parameters1 := [ s = 2 ];
                            [s = 2]

U(c, subs(Parameters1,s));
                              1    
                            - - + 1
                              c    

Your principal question seems to be about the projection. So I'll leave that for "Answers". This Comment is only about your minor issues with labels and captions of the rotated 3D plots.

I just made up this `cast` procedure as a rough way to typeset those fractions in captions. I left out the InvisibleTimes, to pull the 3 and the Pi closer together for the 3*Pi/4.

restart:

cast:=proc(r)
local n,d,t,top;
  n,d:=numer(r),denom(r);
  if type(n,`*`) then
    top:=cat(`mrow(`,cat(`mn(`,convert(op(1,n),string),`)`),
             seq(cat(`,mn(`,convert(op(t,n),string),`)`),t=2..nops(n)));
  else
    top:=cat(`mrow(`,cat(`mn(`,convert(n,string),`)`));
  end if;
  if d<>1 then
    cat(`#mfrac(`,top,`),`,cat(`mn(`,convert(d,string),`)`),`)`);
  else
    cat(`#`,top,`)`,`)`);
  end if;
end proc:

p := plot3d( [ 4+x*cos((1/2)*y), y, x*sin((1/2)*y) ]
  , x = -Pi .. Pi
  , y = 0 .. 2*Pi
  , 'coords' = cylindrical
  , 'style' = patchnogrid
  , 'grid' = [60, 60]
  , 'lightmodel' = light4
  , 'shading' = zhue
  , 'scaling' = constrained
  , 'transparency' = .3
):

P := CodeTools:-Usage( seq(
  plots:-display( plottools:-rotate(p, k*Pi/4, [ [0,0,0], [0,0,1] ]),
                  'axes'='box', 'labels' = [`x`, `y`, `z`],
                  'caption' = [ typeset("3D plot rotated around x-axis by %1",
                                cast(k*Pi/4),
                                         ".") ]
  ) , k = 1 .. 4
) ):

ArraySeqPlot := plots:-display( ArrayTools:-Alias(Array([P]), [round(nops([P])/2), 2],
                                                  C_order) ):

plots:-display(%);

1) It's a hand-wavy way to consider the dy and the dx as separate symbols. In college students are taught that dy/dx is not merely a fraction in the usual sense, but a notation for a concept of differentiation. But then later on they are shown some slick notational abuses/shoftcuts (eg. implicit differentiation) involving "separation" of those symbols. Because it helps solve some examples.

2) Because we're using 1D input mode.

3) There are an untold number of things one can do with maple and a full programming language. The hard-coded palettes and context-menus can only ever be a small portion of that. (Consider, the palettes have no add, just sum, as one trivial case.) you've asked for something nonstandard, so we're pushing the limits.

You wrote that you followed Joe's way, but really you did a mix. You used diff(F(x),x) while Joe deliberately used just the plain ol' fraction dy/dx where the dy and the dx were no special calculus notation (except in our minds' eye).

As shown in my attempt, the `assume positive` can be replaced with `assume x[1]>0 and x[2]>x[1]`. But in order for the RHS (when containing that proper `diff` you and I both used) to integrate Maple likes know that the integrand is continuous. Unfortunately, the 2D input form does not allow the use of extra options to the `int` command. So you might be stuck with the text form int(...,continuous) for the RHS, although the bit inside that could still be in 2D form.

int_v_dsolve-p.mw

First 14 15 16 17 18 19 20 Last Page 16 of 81