acer

32363 Reputation

29 Badges

19 years, 332 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Carl Love Some years ago I wrapped up something very similar into a procedure which I called findroots. I mention it because, for a requested range a..b, it invoked NextZero with an option like maxdistance=b-start rather than just maxdistance=b say.

It looks to me as if your code is intended to restrict further searching if any root is found above bound 10000, because of the while condition on the loop. It doesn't make a huge difference for this example, but it might be that there is some benefit (speed or robustness) by using an option maxdistance= 10000 - R[k-1] rather than maxdistance= 10000 in the NextZero call. The endpoint 10000.0 might need it's own examination, but that's relatively very cheap.

Just an idea.

I think that NextZero has improved over the years. I don't see why fsolve couldn't be amended to call it as a method, possibly while offering the maxsols option for input other than just univariate polynomials. The whole business of repeatedly calling fsolve & its avoid option, as used by Student:-Calculus1:-Roots, is a bit sad.

@Alejandro Jakubi Of course customization is possible, but in my experience there are very many users who want and expect their own favored functionality model to be the default.

When I wrote "thorny dichotomy" I was thinking of the problem that users can and often do want greatly differing and opposing defaults from each other.

@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

First 337 338 339 340 341 342 343 Last Page 339 of 592