acer

32510 Reputation

29 Badges

20 years, 11 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Carl Love I'll note a few things too:

- andmap expects true or false, and will throw an error if the testing procedure such as is returns FAIL before andmap returns. So either a less beautiful test such as t->is(t)=true, or a nested loop, or similar, should be used instead of a raw is call.

- Examples can be found -- including in old posts on this site -- of mixed or unmixed trig, expln, radical, or RootOf expressions whose difference can be reduced to zero using a combination of custom simplification steps, but for which is does not succeed in showing equivalence.

Hence I prefer something more customizable, such as the Normalizer/Testzero method. This can also be set to use `is` of course. It might also be used in general concordance with the pivot checkers, division steps, and row reduction steps done by LUDecomposition for LU and Cholesky methods (because they use it too). I'd favor having the customizability mechanism not be the act of writing new extensions `IsMatrixShape/...`.

I also don't value so much authoring a custom extension to IsMatrixShape, if using it as such an extension requires even more typing. Intead of authoring one's own `IsMatrixShape/IsSymmetric` procedure and having to invoke it by all the keystrokes of IsMatrixShape(A,IsSymmetric) one could more easily author the same procedure and call it IsSymmetric and call it like IsSymmetric(A). I don't see why a sensible and consistent naming convention can't be as useful handling the global namespace as using very long names.

And what does the current implementation of `IsMatrixShape`(M,symmetric) offer that type 'Matrix(symmetric)' does not? If nothing of significance then the former could be sensibly changed to be smarter. I don't see much problem with having it be replaced with something as usefully powerful as a test by normal, which is what Normalizer defaults to. That mentioned typecheck could be used on Matrix A beforehand just as easily as could that call to IsMatrixShape. It's a few keystrokes less to use, since symmetric is not a protected name and hence should be quoted for general safeness.

type(A,'Matrix(symmetric)');
IsMatrixShape(A,'symmetric');

I don't really understand why IsMatrixShape even exists, if its stock methods are not going to be smarter than a corresponding typecheck.

@tomleslie A plain if..then conditional uses evalb as the boolean evaluator. And the help page for evalb documents somewhat well what it means for evalb(A=B) to return true. In the case of general nonnumeric expressions that won't return true unless A and B evaluate to the very same uniquified object in memory. It should therefore not be a surprise that evalb( x*(x-1) = x^2-x ) returns false since those are stored internally by Maple as two distinct symbolic expressions.

On the other hand you have brought up valid concerns. As Carl pointed out, normal may not invoke expand. But it should be recognized that while equivalence testing is a cousin of zero recognition but is not the very same thing. Maple's normal can produce the result that is identically 0 when passed the difference of two mathematically equivalent rational polynomials which might not individually be returned as the very same form by normal. For example,

expr1 := x*(x-1):
expr2 := x^2-x:

normal( expr1 );
                                  x (x - 1)

normal( expr2 );
                                    2    
                                   x  - x

normal( expr2 - expr1 );
                                      0

So throughout your quiz2, say, you'd be on firmer ground by considering expr1-expr2=0 rather than expr1=expr2.

It's also worth considering that people have come down this road many times before. There are several mechanisms for such testing. None is perfect -- not so much because "zero testing is theoretically undecidable", but more practically because the balance between ease-of-us/flexibility and power will never be perfect for all different kinds of Maple user. The following is just an illustration of alternate approaches.

print( eval(Testzero) );
 proc(O) evalb(Normalizer(O) = 0) end proc;

print( eval(Normalizer) );
 proc() option builtin = normal, remember, system;  end proc;

Testzero( expr2 - expr1 );
                                    true

So a potential alternative might be something more configurable. Eg,

restart:

unprotect(`IsMatrixShape/symmetric`);
`IsMatrixShape/symmetric` := proc(M)
   if op([1, 1],M) <> op([1, 2],M) then
       false;
   else;
       andmap(Testzero, M - M^%T );
   end if;
end proc:
protect(`IsMatrixShape/symmetric`);

M1:= Matrix([[0, x*(1+x)],[0, 0]]):
M2:= Matrix([[0, 0],[x+x**2, 0]]):

IsMatrixShape(M1+M2,symmetric);
                              true

Normalizer := t -> simplify(t, trig):

IsMatrixShape(Matrix([[0, sin(x)^2+cos(x)^2],[1, 0]]),symmetric);
                              true

Let's also notice that Maple has a working type-check for symmetric Matrices. So there's not much point in also having a command to do only the same thing, especially if it involves more typing to use. Since the command's name is prefixed by "Is" then why not give it the benefits of the is/assume/assuming mechanism?

# I actually much prefer to have just Testzero be used. I dislike this version
# that I give below because the potentially expensive `is` call cannot be turned off.
# I don't know how to resolve the thorny dichotomy: some users expects correct results
# as often as they can get them, with no surprises. But other others expect a fast and
# efficienct system which can be configured to not waste resources.
restart:

unprotect(`IsMatrixShape/symmetric`);
`IsMatrixShape/symmetric` := proc(M)
   if op([1, 1],M) <> op([1, 2],M) then
       false;
   else;
       andmap(t -> is(Normalizer(t) = 0), M - M^%T );
   end if;
end proc:
protect(`IsMatrixShape/symmetric`);

M1:= Matrix([[0, x*(1+x)],[0, 0]]):
M2:= Matrix([[0, 0],[x+x**2, 0]]):

IsMatrixShape(M1+M2,symmetric);
                              true

IsMatrixShape(Matrix([[0,a],[b,0]]),symmetric) assuming a=b;
                              true

IsMatrixShape(Matrix([[0, sin(x)^2+cos(x)^2],[1, 0]]),symmetric);
                              true

I do believe that merely using EqualEntries is a poor approach, which is what IsMatrixShape(..., symmetric) is doing in Maple 18. It's no better than the simple type-check syntax, I'd say.

restart:
showstat(`IsMatrixShape/symmetric`);

`IsMatrixShape/symmetric` := proc(M)
local S;
   1   if op([1, 1],M) <> op([1, 2],M) then
   2     return false
       end if;
   3   S := Matrix(M,scan = triangular[upper],shape = symmetric);
   4   return EqualEntries(S,M)
end proc

I think it's quite natural to expect at least the following of a symmetric Matrix check:

  • Entrywise comparison as strong as normal by default, and tests the difference against zero, rather than just equality between entries.
  • Configurable, documentable, and consistent with some other parts of Maple.
  • Supports assumptions.

Sure, a new tolerance option for float comparison could also be added to the command. But a configurable mechanism such as provided by Testzero and Normalizer allow that too. They also allow for use of verify,float as well as for more general used of verify for the entrywise scalar comparisons.

Someone might point out that only the upper or lower triangle might need walking, when applying such tests. In the case of using andmap I wonder whether such could be accomplished by that kernel builtin itself.

@Rouben Rostamian  A number of previously valid .mw files would still run with errors in Maple 2015, if the only amendment made to to Maple 2015.0 were to change the behaviour of the F4 action. The issue would still be for backwards compatibility of users' worksheets.

To retain much better compatibility with valid and functioning worksheets produced in M18 and earlier, while still obtaining new convenience (as indicated in the very last example here) it seems to me that the new behaviour could just be amended to:

     A explicit statement terminator such as a colon or semicolon is no longer
     required in 1D Maple Notation input for the last line of an Execution Group.

One way that this popup can occur is if you try and open an empty file (0 bytes) in the Standard GUI, even if it has ".mw" as the file extension.

(And one unfortunate way that zero-length .mw files can be generated is by downloading unzipped .mw attachments from email, using some webmail tools such as client to a misconfigured Outlook server. This because the mail server can misidentify the attachment, since .mw Maple files contain XML.)

So check whether the file is not empty. If it's an .mw file and it's not empty then you might upload it here so that we could try and see what's wrong with it.

acer

Which version of Maple are you using, and which platform (32bit or 64bit?).

acer

@Alejandro Jakubi I think that my narrative here doesn't even come close to WRI's position as described in the link you gave. I'm very much against that view being applied to Maple. And indeed I am one of the relatively few people who has posted fixes to Library code on this site. And many of my Posts here comprise new techniques to make use of Maple's full functionality (whether documented or not). Maple's openness is one of its great strengths.

I merely mentioned waiting for a point-release as something I might do in this particular case. But I also described two approaches, and qualified one of them as being difficult or involved (ie, "tricky"). I haven't advocated that anyone else not try. Everyone is free to decide for themselves how such efforts' cost measure against the benefits.

I see no merit or truth in your Comment. Please snipe in your own Answer, or contribute concretely.

@itsme Well, you might try cycling through all the exports and locals (and similarly for any submodules you discover such as Explore:-Runtime). Fortunately Explore doesn't have a ModuleLoad which redefines any part of itself, so that's one kind of trickiness avoided.

It might be simpler to just try the combination: unprotect, unset kernelopts(opaquemodules) ?, redefine the ModuleApply, reset kernelopts(opaquemodules), and reprotect in an initialization file. And forego resaving to your own .mla archive.

@itsme It is tricky. You need to ensure that the full module, and all its locals, have been read in from system archive. This is not accomplished by just calling eval on the module, and may not be guaranteed just by calling it on a particular example.

Personally, I'd wait for a point-release fix...

@itsme I see. Thanks. I think that it is a bug, and I will submit a bug report.

I agree, the plot array should be aligned (in its sub-table) toward the controls and not away from them. I'd have to dig and debug before figuring out whether a Library-side code change would help -- plot arrays have some funny alignment behavior which is sometimes purely GUI related, and without debugging I'm not sure offhand whether this is a Library or GUI problem.

OTOH, the sub-table that contains the PlotComponent just seems unnecessarily wide. That too might be a fixable bug, possibly by just adjusting some width-related  fudge factor in the code. I'll let you know if I figure something out.

There is also suboptimal behavior in the subtable on the left (where the controller stuff is). The current value (of 1) is in a subtable that is wide enough to contain an "animating" checkbox. But this is not an animated exploration, so the subsubtable that contains the "1" is just unnecessarily wide. I'll submit this too.

Ok, I (sort of) see part of the problem now. The two principal subtables have column weights of 300 and 1210 for the example you gave. That is see by using infolevel[Explore]:=4 I think, It might be possible to fiddle with the weights, but I suppose that changing the alignment of the PlotComponent's subtable might also be ok. It would have to careful, changing the alignment to `left` only when there were no controllers placed on the right (or vice-versa). Perhaps the ideal solution would consist of both kinds of change, relative weights as well as alignment.

BTW, Your example looked "better" in Maple 18.02, but only in the sense that the width of the PlotComponent was being overridden by the outermost Table's behaviour -- fit to GUI 100% of window. This is something related for which I've already logged a Software Change Request -- Explore could get a new width-mode option to control whether the outermost Table was a fixed percentage of GUI window, or a pixel count at default zoom. (The new Tabulate command in Male 2015 has such an option.)

Could you give a (toy, or simple) example where the alignment problem is manifested?

Also, which version of Maple are you using?

acer

@Carl Love Thanks for those timing results. I'd hazard a guess that this (rather high level) parallelism has its performance benefit constrained by shared use of the cache by multiple cores. 

There may be some additional gains possible by some reuse of the column sums in (what we were previously calling `subimage`). If kerneld has horizontal symmetry then the weighted column sums each contribute to a pair of newaa entries, I think. And if the columns of kerneld are scalar multiples of each other then each column sum (scalar) can be reused even more -- it just needs to be adjusted by a scalar weight corresponding to the particular column as it contributes to several newaa[I,j].

OTOH I'm not sure that ImageTools:-Convolution is using the function for that in the Intel IPP library that Maple now uses for SignalProcessing. So may it too could be faster.

Fun stuff.

In Maple 2015.0 the Task multithreaded compiled code I gave immediately above runs about twice as fast as in Maple 18.02 on my 64bit Windows 7, for that 600x800 example.

@Carl Love You should check your results for correctness, when using `Compile` on a procedure containing calls to `add` over float values. It's a "known bug" that it often get confused about the types in play, in such a mix. As you originally have it the code intended to be compiled will not produce correct results, I suspect. The workaround of replacing `add` calls with explicit loops may slow down under evalhf. So one may need distinct versions for running compiled versus running under evalhf.

The following uses explicit loops instead of `add`, for the Compile version. The temporary accumulator for the sum is a scalar variable, since using `newaa[i,j]` instead as accumulator can (MS-Windows, LLVM compiler) incur a big performance hit due to unnecessary referencing into array memory.

Also, as mentioned before, I also obtain another performance gain by computing i-k-1+m outside of the n-loop (stored in `mindex`).

And lastly, by having the procedure `Convolve_C` accept the lower and upper index values for the i-loop as arguments, the Threads:-Task model can be applied to parallelize the computation.

 

restart:

kernelopts(version);
kernelopts(numcpus);

          Maple 18.02, X86 64 WINDOWS, Oct 20 2014, Build ID 991181
                                      8

k := 3:

filelocation := cat(kernelopts(datadir),"/ship.jpg"):
zimage := ImageTools:-Read(filelocation):
zwidth := ImageTools:-Width(zimage):
zHeight := ImageTools:-Height(zimage):
#zimage := ImageTools:-Scale(zimage, 1..2*zHeight, 1..2*zwidth):

kernell := 2*k+1:
kerneld := Matrix(kernell, kernell, fill=1/kernell^2,
                  datatype=float[8], order=C_order):

st := time[real]():

new1zumage := ImageTools:-Convolution(zimage, kerneld):

Time[Convolution] := time[real]()-st;

                         Time[Convolution] := 0.112

aa := zimage(1..,1..,1):
bb := zimage(1..,1..,2):
cc := zimage(1..,1..,3):

newaa := copy(aa):
newaa1 := copy(aa):
newaa2 := copy(aa):
newaa3 := copy(aa):

R := LinearAlgebra:-RowDimension(aa);
C := LinearAlgebra:-ColumnDimension(aa);

                                  R := 400
                                  C := 600

st := time[real]():
for i from k+1 to R-k-1 do
     for j from k+1 to C-k-1 do
          newaa[i,j] := add(x, x = aa(i-k..i+k, j-k..j+k) *~ kerneld);
     end do;
end do;
Time[NativeMaple] := time[real]()-st;

                         Time[NativeMaple] := 16.392

Convolve := proc(aa::Array(datatype=float[8],order=C_order),

                 newaa::Array(datatype=float[8],order=C_order),

                 kerneld::Array(datatype=float[8],order=C_order),

                 k::posint,R::posint,C::posint)
  local i, j, m, n, mindex,
        `k+1` := k+1, `C-k-1` := C-`k+1`, `R-k-1` := R-`k+1`,
        kernell := 2*k+1, `i-k-1`, `j-k-1`;
  for i from `k+1` to `R-k-1` do
    `i-k-1` := i-`k+1`;
    for j from `k+1` to `C-k-1` do
      `j-k-1` := j-`k+1`;
      newaa[i,j] := add(add(aa[`i-k-1`+m,`j-k-1`+n]*kerneld[m,n],
                            n=1..kernell),m=1..kernell);
    end do;
  end do;
end proc:

st := time[real]():
evalhf(Convolve(aa, newaa1, kerneld, k, R, C)):
Time[evalhf] := time[real]()-st;
LinearAlgebra:-Norm( Matrix(newaa)-Matrix(newaa1) );

                            Time[evalhf] := 7.612
                                                -14
                          9.69779812010074238 10   

ConvolveC := proc(aa::Array(datatype=float[8],order=C_order),
                  newaa::Array(datatype=float[8],order=C_order),
                  kerneld::Array(datatype=float[8],order=C_order),
                  k::posint,R::posint,C::posint,
                  rlo::posint,rhi::posint)
  option threadsafe;
  local i::integer[4], j::integer[4], m::integer[4], n::integer[4],
        mindex::integer[4], kp1::integer[4], Cmkp1::integer[4],
        kernell::integer[4], imkp1::integer[4], jmkp1::integer[4],
        temp1::float[8];
  kp1 := k+1;
  kernell := k+kp1;
  Cmkp1 := C-kp1;
  for i from rlo to rhi do
    imkp1 := i-kp1;
    for j from kp1 to Cmkp1 do
      jmkp1 := j-kp1;
      temp1 := 0.0;
      for m from 1 to kernell do
        mindex := imkp1+m;
        for n from 1 to kernell do
          temp1 := temp1 + aa[mindex,jmkp1+n] * kerneld[m,n];
        end do;
      end do;
      newaa[i,j] := temp1;
    end do;
  end do;
end proc:

Convolve_C := Compiler:-Compile(ConvolveC):

st := time[real]():
Convolve_C(aa, newaa2, kerneld, k, R, C, k+1, R-(k+1)):
Time[compiled] := time[real]()-st;
LinearAlgebra:-Norm( Matrix(newaa)-Matrix(newaa2) );

                           Time[compiled] := 0.509

                                                -14
                          8.53206394424432801 10   

Convolve_Task := proc(aa, newaa, kerneld, k, R, C, rlo, rhi)
local rmid;
  if 11 < rhi - rlo then
    rmid := floor(1/2*rhi-1/2*rlo) + rlo;
    Threads:-Task:-Continue(':-null',
                            ':-Task'=[Convolve_Task,aa,newaa,kerneld,k,R,C,rlo,rmid],
                            ':-Task'=[Convolve_Task,aa,newaa,kerneld,k,R,C,rmid+1,rhi])
  else
    Convolve_C(aa,newaa,kerneld,k,R,C,rlo,rhi);
  end if;
end proc:

st := time[real]():
Threads:-Task:-Start(Convolve_Task, aa, newaa3, kerneld, k, R, C, k+1, R-(k+1));
Time[Compiled_parallel] := time[real]()-st;
LinearAlgebra:-Norm( Matrix(newaa)-Matrix(newaa3) );

                      Time[Compiled_parallel] := 0.232
                                                -14
                          8.53206394424432801 10   

Times = Vector(sort(op(eval(Time)),(a,b)->rhs(a)>rhs(b)));

                             [  NativeMaple = 16.392   ]
                             [                         ]
                             [     evalhf = 7.612      ]
                             [                         ]
                     Times = [    compiled = 0.509     ]
                             [                         ]
                             [Compiled_parallel = 0.232]
                             [                         ]
                             [   Convolution = 0.112   ]

 

 

Download convolve_compiled3.mw

@Carl Love Yes, that top-level code isn't the same as calling `Convolve` without evalhf. Rather, it is an elementwise (vectorized) use of kernel builtin `*` on float[8] rtables which it knows how to treat well. So, yes, it's already parallelized in another way, and is fast despite the fact that it repeatedly creates `subimage` as a collectible garbage rtable.

Attached is a sheet with a few ideas on speedup. This does not include what I'd expect to be best: Compiled and multithreaded with the Task mode.

The top level code can be optimized at its second stage, by rolling the two nested `add` calls into just one call.

Yes, I believe that using order=C_order on `subimage` and `kerneld` can make it a bit faster -- as long as loops or the nested `add` calls walk along rows at their innermost layer (so as to walk adjacent entries in memory).

More saving can be had by constucting part of the `subimage` indicies, for re-use, at the start of each nested i and j loop. (See `ipart`, `jpart`, and `mindex`.)

And then one can notice that `subimage` is not really needed in such a proc.

restart:

kernelopts(version);

          Maple 18.02, X86 64 WINDOWS, Oct 20 2014, Build ID 991181

k := 3:

filelocation := cat(kernelopts(datadir),"/ship.jpg"):

zimage := ImageTools:-Read(filelocation):

zwidth := ImageTools:-Width(zimage):

zHeight := ImageTools:-Height(zimage):
zimage := ImageTools:-Scale(zimage, 1..2*zHeight, 1..2*zwidth):

kernell := 2*k+1:

kerneld := Matrix(kernell, kernell, fill=1/kernell^2, datatype=float[8], order=C_order):

st:= time():

new1zumage:= ImageTools:-Convolution(zimage,kerneld):

time()-st;

                                    0.203

#ImageTools:-View(new1zumage);

aa:= zimage(1..,1..,1):

bb:= zimage(1..,1..,2):

cc:= zimage(1..,1..,3):

subimage := Matrix(kernell, kernell, datatype=float[8], order=C_order):

newaa := copy(aa):
newaa1 := copy(aa):
newaa2 := copy(aa):
newaa3 := copy(aa):

R := LinearAlgebra:-RowDimension(aa);

C := LinearAlgebra:-ColumnDimension(aa);

                                  R := 800
                                  C := 1200

st := time[real]():

for i from k+1 to R-k-1 do

     for j from k+1 to C-k-1 do

          subimage := aa(i-k..i+k, j-k..j+k) *~ kerneld;
          newaa[i,j] := add(s, s=subimage);

     end do

end do;

T1s := time[real]()-st;

                                T1s := 40.814

st := time[real]():

for i from k+1 to R-k-1 do

     for j from k+1 to C-k-1 do

          subimage := aa(i-k..i+k,j-k..j+k) *~ kerneld;

          newaa[i,j] := add(add(subimage[m,n], m=1..kernell), n=1..kernell)

     end do

end do;

T1:= time[real]()-st;

                                T1 := 54.961

Convolve:= proc(aa, newaa, kerneld, subimage, k, R, C)

local i, j, m, n, `k+1` := k+1, `C-k-1` := C-`k+1`, kernell := 2*k+1;     

     for i from `k+1` to R-k-1 do

          for j from `k+1` to `C-k-1` do

               for m to kernell do

                    for n to kernell do

                         subimage[m,n] := aa[i-k+m-1, j-k+n-1] * kerneld[m,n]

                    end do

               end do;

               newaa[i,j] := add(add(subimage[m,n], m=1..kernell), n=1..kernell)

          end do

     end do

end proc:

st := time[real]():

evalhf(Convolve(aa, newaa1, kerneld, subimage, k, R, C)):

T2 := time[real]()-st;

                                T2 := 28.373

Convolve2:= proc(aa, newaa, kerneld, subimage, k, R, C)

local i, j, m, n, `k+1`:= k+1, `C-k-1` := C-`k+1`, kernell := 2*k+1,
      ipart, jpart, mindex;     

     for i from `k+1` to R-k-1 do
          ipart := i-k-1;

          for j from `k+1` to `C-k-1` do
               jpart := j-k-1;

               for m to kernell do
                    mindex := ipart+m;

                    for n to kernell do

                         subimage[m,n] := aa[mindex, jpart+n] * kerneld[m,n]

                    end do

               end do;

               newaa[i,j] := add(add(subimage[m,n], n=1..kernell), m=1..kernell)

          end do

     end do

end proc:

st := time[real]():

evalhf(Convolve2(aa, newaa2, kerneld, subimage, k, R, C)):

T3 := time[real]()-st;

                                T3 := 24.766

Convolve3:= proc(aa, newaa, kerneld, k, R, C)

local i, j, m, n, `k+1`:= k+1, `C-k-1` := C-`k+1`, kernell := 2*k+1,
      ipart, jpart;     

      for i from `k+1` to R-k-1 do
          ipart := i-k-1;

          for j from `k+1` to `C-k-1` do
               jpart := j-k-1;

               newaa[i,j] := add(add(aa[ipart+m, jpart+n] * kerneld[m,n],
                                     n=1..kernell), m=1..kernell)

          end do

     end do

end proc:

st := time[real]():

evalhf(Convolve3(aa, newaa3, kerneld, k, R, C)):

T4 := time[real]()-st;

                                T4 := 19.051

LinearAlgebra:-Norm(Matrix(newaa)-Matrix(newaa1));
LinearAlgebra:-Norm(Matrix(newaa)-Matrix(newaa2));
LinearAlgebra:-Norm(Matrix(newaa)-Matrix(newaa3));

                                     0.
                                                -14
                          6.30606677987088915 10   
                                                -14
                          6.30606677987088915 10   

 


Download convolve_evalhf.mw

@rollermonkey Your claims about my suggestions appear to be entirely false, and appear to depend on your continued mistyping of key syntax, or misspelling (you now appear to have done with(Plots) not with(plots) as I suggested), or perhaps in part due to a faulty installation. Carl's Answer covers it in detail, I see.

First 340 341 342 343 344 345 346 Last Page 342 of 595