Carl Love

Carl Love

16959 Reputation

24 Badges

7 years, 234 days
Mt Laurel, New Jersey, United States
My name was formerly Carl Devore. I was in the PhD math program at University of Delaware until 2005. I was very active in the Maple community at that time.

MaplePrimes Activity


These are replies submitted by Carl Love

@Carl Love I pared my code down to a single regular-expression substitution:

showstat_no_nums:= p->
    `debugger/printf`(
        "%s",  
        StringTools:-RegSubs(
            "\n[ 1-9][^\n][^\n][^\n]"= "\n ", 
            debugopts('procpretty'= p)
        )
    )
:

 

@emendes 

Sorry for the delay in answering. It was because both of your questions require lengthy and technical answers.
 

Re: "Memory used" as reported by CodeTools:-Usage: This number is reported with a dimension of information (bytes or multiples thereof). This is incorrect! What this number actually represents is Int(M(t), t= time program runs) where M(t) is the amount of memory "being actively used" at time t. Thus the correct dimensions are information*time (e.g., GiB-s (gibibyte-seconds)). It is sad and somewhat shocking that no-one at MapleSoft has noticed and corrected this. (The same dimension error is incorrectly documented for kernelopts(bytesused).) Many systems and computer scientists avoid the confusion about this issue by naming this measurement "memory integral" rather than "memory used" and by using the correct dimensions.

I don't know the precise definition of "being actively used"; I am trying to figure that out from experimentation. There is also the issue of how time is measured. Is it real time? cpu time? something else? This is a bit easier to determine by experiment, and my best guess so far is that it's pretty close to kernelopts(cputime) - kernelopts(gctotaltime).

You're justifiably worried about running out of memory when doing the larger tuples. The number that you need to watch is not reported by CodeTools:-Usage. It does however appear on the worksheet status line (the bottom right), along with cputime.

Note that memory is nearly constantly being recycled by garbage collection, so the same chunk of memory may be "used" many times during a program run. Each such re-use adds to the incorrectly named "memory used".
 

Re: option remember and option cache: Both of these are attempts to save the cputime needed to re-call a procedure with arguments that it's already seen at the expense of the memory needed to store those results. The process is technically called memoization in computer science literature.

Consider the subprocedure of OrbitPartition. It may be called millions of times during an execution of the program, but it only has exactly 13 possible pairs of arguments. Thus the memory needed to store the pre-computed results is trivial. Using option remember means any previously computed results will be remembered for as long as the procedure is in the visible scope.

Now consider the procedure ClassifyI in the setup code. It also may be called millions or times, but with millions of different arguments. But it's called twice in succession with the same argument---once in AllowedRep and once in DisallowedOrbit. By using option cache, a scrolling cache of the most recently used results is maintained. You can specify the maximum number of entries in this cache in the option. In this case, option cache(1) would work. The only reason that I didn't put a number is that a default-size cache might supply some useful debugging information if that were needed. The cache or remember table of a procedure P (it can't have both, but it may have neither) can be directly inspected via op(4, eval(P)).
 

Re: Making it a module: Do you want the setup specifications inside (local to) the module?

@JAMET Sorry about that. I had just noticed that myself, and I uploaded a new worksheet a few minutes ago. Try it again.

I have often had trouble with Grid:-Seq when the chunks were too small. I'd recommend that each chunk take at least 2 seconds to process. Your 500 seems both rather arbitrary and way too small.

I don't think that it can be done unless you know the positions of the vertices. If you do know (for example, if you can measure them from an existing picture), then you can use SetVertxPositions. I doubt that there is any algorithm to do it without knowing those vertex positions.

To give you an idea of the problems with constructing an algorithm, consider the triangle inequality.

@mmcdara The output of showstat is equivalent to the output of printf; so, yes, it's "side-effect" output rather than "return value" output. That's why I went one step back through the showstat code to get the output from debugopts(procpretty= ...). This is identical to the standard showstat output, but in string (actually name) form.

@Pascal4QM Thanks. As you can probably tell, I worked it out by deconstructing showstat itself.

@acer You can use a line-width specifier in the format code such as "%40P". Then, the lines will be sensibly broken and indented with a maximum length of 40. It's quite impressive how nicely it's done.

@dharr 

dsolve(..., numeric) does try to solve ODEs for their highest derivatives before proceeding. But if there are multiple solutions, even trivial ones, it refuses to arbitrarily select one.

There are other issues that may hamper the attempt that are relatively easy to address "by hand". For example, in the case at hand, I needed to differentiate one of the equations. That often necessitates additional initial\boundary conditions, but it didn't in this case.

The physics of the original equations is wrong. The algebra by which I first used solve to get dsolve to solve those equations is very standard. Other than correcting g to negative, I didn't look at the physics of the equations; I figured that that'd be a job better suited for you. Indeed, I almost concluded the Answer with "I'll let Rouben address the physics."

@acer The OP is attemptimg to solve this problem, which was posted by another user a few days ago. I guess that they're taking the same course. Thus, the exponent in the integrand is intended to be 3/2, not 3. Maple can still solve it symbolically, although it's not clear whether that's desirable.

@emendes 

Here is my code that processes all 4 extra conditions. I get 123 3-tuples, just as you did. 

restart:

gc()
;

OrbitPartition:= (
    parms::And(set, listlist &under [op]),
    tuple_size::posint,
    permutations_by_index::list(table),
    AllowedRep::appliable,    
    DisallowedOrbit::appliable
)->
local
    C, i, Done:= table(),
    T:= proc(T,j) option remember; `if`(assigned(T[j]), T[j], j) end proc,
    conds:= varCoef-> map(x-> T~(permutations_by_index, x), varCoef),
    Orbit:= x->
        local
            r:= x,
            R:= {
                do  
                    Done[r]:=
                        #Caution: AllowedRep must be true in order to call
                        #    DisallowedOrbit without error.
                        if AllowedRep(r) then if DisallowedOrbit(r) then return else r fi
                        else
                        fi
                until assigned(Done[(r:= conds(r))])
            }
        ;
        if R <> {} then R[1] fi
;   
    {seq}(
        if assigned(Done[(i:= parms[[seq(C+~1)]])]) then else Orbit(i) fi,
        C= Iterator:-Combination(nops(parms), tuple_size)
    )
:

#Problem setup:
#--------------
#
nIs:= 3: Is:= {$1..nIs}:
parms:= {seq(seq([i,j], j= 0..9), i= Is)}:
TJ:= table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]):
TI:= table([2= 3, 3= 2]):
#
#Exclusion conditions:
#---------------------
#
ClassifyI:= proc(S) option cache; ListTools:-Classify(index, S, 1) end proc
:
AllowedRep:= S->
        Is = op~(1,S)                                  #condition 1
    and max(op~(2,S)) >= 4                             #condition 2
    and nops({entries((op~)~(2,ClassifyI(S)))}) = nIs  #condition 3
:
DisallowedOrbit:= S->
local C:= ClassifyI(S);
    #Caution: Condition 1 must be true in order to check condition 4 without error.
    #condition 4:
    ormap(k-> op~(2,C[k]) subset [{0,1,4}, {0,2,7}, {0,3,9}][k], Is)
:   

#Test run:
newabc:= CodeTools:-Usage(
    OrbitPartition(
        parms, 3, [TI,TJ],
        #(*        
        AllowedRep, DisallowedOrbit
        #*)
        (*
        ()-> true, ()-> false
        *)
    )
        
):

memory used=22.23MiB, alloc change=0 bytes, cpu time=172.00ms, real time=141.00ms, gc time=62.50ms

n:= nops(newabc);

123

interface(rtablesize= (n4:= iquo(n,4))):

map2(`?[]`, <(`<|>`@op)~([newabc[]])[]>, [[k*n4+1..(k+1)*n4] $ k=0..3, [4*n4+1..-1]])[];

Matrix(30, 3, {(1, 1) = [1, 2], (1, 2) = [2, 1], (1, 3) = [3, 4], (2, 1) = [1, 2], (2, 2) = [2, 1], (2, 3) = [3, 5], (3, 1) = [1, 2], (3, 2) = [2, 1], (3, 3) = [3, 6], (4, 1) = [1, 2], (4, 2) = [2, 1], (4, 3) = [3, 7], (5, 1) = [1, 2], (5, 2) = [2, 1], (5, 3) = [3, 8], (6, 1) = [1, 2], (6, 2) = [2, 3], (6, 3) = [3, 4], (7, 1) = [1, 2], (7, 2) = [2, 3], (7, 3) = [3, 5], (8, 1) = [1, 2], (8, 2) = [2, 3], (8, 3) = [3, 6], (9, 1) = [1, 2], (9, 2) = [2, 3], (9, 3) = [3, 7], (10, 1) = [1, 2], (10, 2) = [2, 3], (10, 3) = [3, 8], (11, 1) = [1, 2], (11, 2) = [2, 4], (11, 3) = [3, 1], (12, 1) = [1, 2], (12, 2) = [2, 4], (12, 3) = [3, 5], (13, 1) = [1, 2], (13, 2) = [2, 4], (13, 3) = [3, 6], (14, 1) = [1, 2], (14, 2) = [2, 4], (14, 3) = [3, 7], (15, 1) = [1, 2], (15, 2) = [2, 4], (15, 3) = [3, 8], (16, 1) = [1, 2], (16, 2) = [2, 5], (16, 3) = [3, 1], (17, 1) = [1, 2], (17, 2) = [2, 5], (17, 3) = [3, 4], (18, 1) = [1, 2], (18, 2) = [2, 5], (18, 3) = [3, 6], (19, 1) = [1, 2], (19, 2) = [2, 5], (19, 3) = [3, 7], (20, 1) = [1, 2], (20, 2) = [2, 5], (20, 3) = [3, 8], (21, 1) = [1, 2], (21, 2) = [2, 6], (21, 3) = [3, 1], (22, 1) = [1, 2], (22, 2) = [2, 6], (22, 3) = [3, 4], (23, 1) = [1, 2], (23, 2) = [2, 6], (23, 3) = [3, 5], (24, 1) = [1, 2], (24, 2) = [2, 6], (24, 3) = [3, 7], (25, 1) = [1, 2], (25, 2) = [2, 6], (25, 3) = [3, 8], (26, 1) = [1, 2], (26, 2) = [2, 8], (26, 3) = [3, 1], (27, 1) = [1, 2], (27, 2) = [2, 8], (27, 3) = [3, 4], (28, 1) = [1, 2], (28, 2) = [2, 8], (28, 3) = [3, 5], (29, 1) = [1, 2], (29, 2) = [2, 8], (29, 3) = [3, 6], (30, 1) = [1, 2], (30, 2) = [2, 8], (30, 3) = [3, 7]}), Matrix(30, 3, {(1, 1) = [1, 2], (1, 2) = [2, 9], (1, 3) = [3, 1], (2, 1) = [1, 2], (2, 2) = [2, 9], (2, 3) = [3, 4], (3, 1) = [1, 2], (3, 2) = [2, 9], (3, 3) = [3, 5], (4, 1) = [1, 2], (4, 2) = [2, 9], (4, 3) = [3, 6], (5, 1) = [1, 2], (5, 2) = [2, 9], (5, 3) = [3, 7], (6, 1) = [1, 2], (6, 2) = [2, 9], (6, 3) = [3, 8], (7, 1) = [1, 5], (7, 2) = [2, 1], (7, 3) = [3, 2], (8, 1) = [1, 5], (8, 2) = [2, 1], (8, 3) = [3, 4], (9, 1) = [1, 5], (9, 2) = [2, 1], (9, 3) = [3, 6], (10, 1) = [1, 5], (10, 2) = [2, 1], (10, 3) = [3, 7], (11, 1) = [1, 5], (11, 2) = [2, 1], (11, 3) = [3, 8], (12, 1) = [1, 5], (12, 2) = [2, 3], (12, 3) = [3, 1], (13, 1) = [1, 5], (13, 2) = [2, 3], (13, 3) = [3, 2], (14, 1) = [1, 5], (14, 2) = [2, 3], (14, 3) = [3, 4], (15, 1) = [1, 5], (15, 2) = [2, 3], (15, 3) = [3, 6], (16, 1) = [1, 5], (16, 2) = [2, 3], (16, 3) = [3, 7], (17, 1) = [1, 5], (17, 2) = [2, 3], (17, 3) = [3, 8], (18, 1) = [1, 5], (18, 2) = [2, 4], (18, 3) = [3, 1], (19, 1) = [1, 5], (19, 2) = [2, 4], (19, 3) = [3, 2], (20, 1) = [1, 5], (20, 2) = [2, 4], (20, 3) = [3, 6], (21, 1) = [1, 5], (21, 2) = [2, 4], (21, 3) = [3, 7], (22, 1) = [1, 5], (22, 2) = [2, 4], (22, 3) = [3, 8], (23, 1) = [1, 5], (23, 2) = [2, 6], (23, 3) = [3, 1], (24, 1) = [1, 5], (24, 2) = [2, 6], (24, 3) = [3, 2], (25, 1) = [1, 5], (25, 2) = [2, 6], (25, 3) = [3, 4], (26, 1) = [1, 5], (26, 2) = [2, 6], (26, 3) = [3, 7], (27, 1) = [1, 5], (27, 2) = [2, 6], (27, 3) = [3, 8], (28, 1) = [1, 5], (28, 2) = [2, 8], (28, 3) = [3, 1], (29, 1) = [1, 5], (29, 2) = [2, 8], (29, 3) = [3, 2], (30, 1) = [1, 5], (30, 2) = [2, 8], (30, 3) = [3, 4]}), Matrix(30, 3, {(1, 1) = [1, 5], (1, 2) = [2, 8], (1, 3) = [3, 6], (2, 1) = [1, 5], (2, 2) = [2, 8], (2, 3) = [3, 7], (3, 1) = [1, 5], (3, 2) = [2, 9], (3, 3) = [3, 1], (4, 1) = [1, 5], (4, 2) = [2, 9], (4, 3) = [3, 2], (5, 1) = [1, 5], (5, 2) = [2, 9], (5, 3) = [3, 4], (6, 1) = [1, 5], (6, 2) = [2, 9], (6, 3) = [3, 6], (7, 1) = [1, 5], (7, 2) = [2, 9], (7, 3) = [3, 7], (8, 1) = [1, 5], (8, 2) = [2, 9], (8, 3) = [3, 8], (9, 1) = [1, 7], (9, 2) = [2, 1], (9, 3) = [3, 2], (10, 1) = [1, 7], (10, 2) = [2, 1], (10, 3) = [3, 4], (11, 1) = [1, 7], (11, 2) = [2, 1], (11, 3) = [3, 5], (12, 1) = [1, 7], (12, 2) = [2, 1], (12, 3) = [3, 6], (13, 1) = [1, 7], (13, 2) = [2, 1], (13, 3) = [3, 8], (14, 1) = [1, 7], (14, 2) = [2, 3], (14, 3) = [3, 1], (15, 1) = [1, 7], (15, 2) = [2, 3], (15, 3) = [3, 2], (16, 1) = [1, 7], (16, 2) = [2, 3], (16, 3) = [3, 4], (17, 1) = [1, 7], (17, 2) = [2, 3], (17, 3) = [3, 5], (18, 1) = [1, 7], (18, 2) = [2, 3], (18, 3) = [3, 6], (19, 1) = [1, 7], (19, 2) = [2, 3], (19, 3) = [3, 8], (20, 1) = [1, 7], (20, 2) = [2, 4], (20, 3) = [3, 1], (21, 1) = [1, 7], (21, 2) = [2, 4], (21, 3) = [3, 2], (22, 1) = [1, 7], (22, 2) = [2, 4], (22, 3) = [3, 5], (23, 1) = [1, 7], (23, 2) = [2, 4], (23, 3) = [3, 6], (24, 1) = [1, 7], (24, 2) = [2, 4], (24, 3) = [3, 8], (25, 1) = [1, 7], (25, 2) = [2, 5], (25, 3) = [3, 1], (26, 1) = [1, 7], (26, 2) = [2, 5], (26, 3) = [3, 2], (27, 1) = [1, 7], (27, 2) = [2, 5], (27, 3) = [3, 4], (28, 1) = [1, 7], (28, 2) = [2, 5], (28, 3) = [3, 6], (29, 1) = [1, 7], (29, 2) = [2, 5], (29, 3) = [3, 8], (30, 1) = [1, 7], (30, 2) = [2, 6], (30, 3) = [3, 1]}), Matrix(30, 3, {(1, 1) = [1, 7], (1, 2) = [2, 6], (1, 3) = [3, 2], (2, 1) = [1, 7], (2, 2) = [2, 6], (2, 3) = [3, 4], (3, 1) = [1, 7], (3, 2) = [2, 6], (3, 3) = [3, 5], (4, 1) = [1, 7], (4, 2) = [2, 6], (4, 3) = [3, 8], (5, 1) = [1, 7], (5, 2) = [2, 8], (5, 3) = [3, 1], (6, 1) = [1, 7], (6, 2) = [2, 8], (6, 3) = [3, 2], (7, 1) = [1, 7], (7, 2) = [2, 8], (7, 3) = [3, 4], (8, 1) = [1, 7], (8, 2) = [2, 8], (8, 3) = [3, 5], (9, 1) = [1, 7], (9, 2) = [2, 8], (9, 3) = [3, 6], (10, 1) = [1, 7], (10, 2) = [2, 9], (10, 3) = [3, 1], (11, 1) = [1, 7], (11, 2) = [2, 9], (11, 3) = [3, 2], (12, 1) = [1, 7], (12, 2) = [2, 9], (12, 3) = [3, 4], (13, 1) = [1, 7], (13, 2) = [2, 9], (13, 3) = [3, 5], (14, 1) = [1, 7], (14, 2) = [2, 9], (14, 3) = [3, 6], (15, 1) = [1, 7], (15, 2) = [2, 9], (15, 3) = [3, 8], (16, 1) = [1, 8], (16, 2) = [2, 1], (16, 3) = [3, 2], (17, 1) = [1, 8], (17, 2) = [2, 1], (17, 3) = [3, 4], (18, 1) = [1, 8], (18, 2) = [2, 1], (18, 3) = [3, 5], (19, 1) = [1, 8], (19, 2) = [2, 1], (19, 3) = [3, 6], (20, 1) = [1, 8], (20, 2) = [2, 1], (20, 3) = [3, 7], (21, 1) = [1, 8], (21, 2) = [2, 3], (21, 3) = [3, 2], (22, 1) = [1, 8], (22, 2) = [2, 3], (22, 3) = [3, 4], (23, 1) = [1, 8], (23, 2) = [2, 3], (23, 3) = [3, 5], (24, 1) = [1, 8], (24, 2) = [2, 3], (24, 3) = [3, 6], (25, 1) = [1, 8], (25, 2) = [2, 3], (25, 3) = [3, 7], (26, 1) = [1, 8], (26, 2) = [2, 4], (26, 3) = [3, 5], (27, 1) = [1, 8], (27, 2) = [2, 4], (27, 3) = [3, 6], (28, 1) = [1, 8], (28, 2) = [2, 4], (28, 3) = [3, 7], (29, 1) = [1, 8], (29, 2) = [2, 5], (29, 3) = [3, 6], (30, 1) = [1, 8], (30, 2) = [2, 5], (30, 3) = [3, 7]}), Matrix(3, 3, {(1, 1) = [1, 8], (1, 2) = [2, 6], (1, 3) = [3, 5], (2, 1) = [1, 8], (2, 2) = [2, 6], (2, 3) = [3, 7], (3, 1) = [1, 8], (3, 2) = [2, 9], (3, 3) = [3, 7]})

 

Download Symmetries.mw

Answers (in no particular order):

  1. Re: Why op rather than `?[]`: The syntax `?[]`(A, [k]) requires [k], whereas op(k, A) does not. Those extra [ ] contribute to the garbage that must be periodically collected. Garbage collection is a major user of processor time.
  2. Re: Condition 3 and tuple length: I misunderstood condition 3. I understand it now, and you're correct, the upper limit of tuple length is much greater than the number of distinct js. I needed to see both of your examples to understand condition 3.
  3. Re: Maximum orbit length: Since conds(conds(x)) = x for any (using the current T1 and T2), the maximum orbit length is 2. By the definition of orbit, the orbit length is the number of applications of conds needed to get back to the original. Suppose that T1 were changed to table([1= 2, 2= 3, 3= 1]). Then there would be length-3 orbits. Can this ever happen? This is equivalent to asking Can a variable ever have more than one "counterpart"?

@AHSAN Here's a Maple worksheet: NewtonRaphson.mw

@emendes 

Questions: 

  1. For your new (4th) condition, you said "..., as well as its counterpart". Does "..., as well as its counterpart" also apply to the other three conditions? In other words, is it necessarily the case for all possible conditions that the invalidity of one member of its orbit makes the whole orbit invalid even though some members of the orbit may not violate the condition? If so, I'll need to change the above code.
  2. Using tables T1 and T2 as above, conds is an involution, i.e., all nontrivial orbits are length 2. Will that always be true, or is that just a property of this particular T1 and T2? If it's always true, I may be able to make the code more efficient.

Answers:

  1. The 3rd condition: You show an example of failure of the 3rd condition (i.e., an invalid set), not satisfaction (i.e., a valid set). Yes, of course, invalid sets can have an arbitrarily high number of elements; however, valid sets necessarily have a number of elements at most as large as the number of distinct js. In our case, that means that it's impossible for an 11-or-greater-tuple to be valid.
  2. `?[]` vs op: The function `?[]`(A, [k]) is the prefix functional form of A[k] (in all cases). For our cases, it can be (and should be!) replaced by op(k, A). It's my fault for not immediately realizing that on my own. My next version will use op instead of `?[]` wherever possible.

@emendes You can check for those additional conditions like this (this is coded for Maple 2019 or later):

OrbitPartition:= (
    parms::And(set, listlist &under [op]),
    tuple_size::posint,
    permutations_by_index::list(table),
    Allowed::appliable
)->
local 
    C, i, Done:= table(),
    T:= proc(T,j) option remember; `if`(assigned(T[j]), T[j], j) end proc,
    conds:= varCoef-> map(x-> T~(permutations_by_index, x), varCoef),
    Orbit:= x-> 
        local 
            r:= x, 
            R:= {
                do 
                    Done[r]:= if Allowed(r, tuple_size) then r else fi 
                until assigned(Done[(r:= conds(r))])
            }
        ;
        if R <> {} then R[1] fi
;   
    {seq}(
        if assigned(Done[(i:= parms[[seq(C+~1)]])]) then else Orbit(i) fi,
        C= Iterator:-Combination(nops(parms), tuple_size)
    )
:
parms:= {seq(seq([i,j], j= 0..9), i= 1..3)}:
T1:= table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]):
T2:= table([2= 3, 3= 2]):
Allowed:= (S,sz)->
    local J:= map(`?[]`, S, [2]); 
    {1,2,3} subset map(`?[]`, S, [1]) and max(J) >= 4 and nops(J)=sz
:
   
newabc:= CodeTools:-Usage(OrbitPartition(parms, 3, [T2,T1], Allowed)): 
memory used=10.85MiB, alloc change=-9.46MiB, cpu time=235.00ms, 
real time=137.00ms, gc time=156.25ms

nops(newabc);
                              358

Your 3rd condition cannot possibly be satisfied by any set with more than 10 elements, so there's no need to check those tuple sizes.

5 6 7 8 9 10 11 Last Page 7 of 513