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

You want to "simplify" the given irrational algebraic constant? It's irrational, so it doesn't get much simpler than it already is. It's equal to 5*2^(-59/6), but that's hardly any simpler. What the point? If you're solving a practical physics problem about dimensioned quantities (such as temperature), then exact answers (other than 0, infinity, etc.) are irrelevant. Just use evalf.

@vv Yes, it's an impressive program. I guess that the input inequalities must be converted to polynomial? If that's the case, you may be able to avoid most solve failures by using solve(..., parametric, real).

@JAMET Now verifying the identity is trivial by hand or with Maple:

is((z1+z2)^2/(z1*z2) = (z1/z2 + z2/z1 + 2));
             
true

 

@JAMET Let z1 and z2 both equal 1. Then the left side is 2 and the right side is 4. They're not equal.

@JAMET I think that you transcribed something wrong: (z1+z2)^2/(z1+z2) = (z1+z2), obviously.

@emendes Here is new code that processes all 6 conditions. You'll see that I added a mechanism to make it easy to exclude conditions (3 is excluded here) and to specify their order (I used 1,2,5,4,6).

I also included to kludge to workaround the error that occurs immediately after restart.

OrbitPartition:= module()
option
    `Conceptual author: emendes`,
    `Maple code author: Carl Love <carl.j.love@gmail.com> 2020-May-7`
;
uses It= Iterator;
local 
    #Problem setup:
    #==============
    nIs:= 3, #number of distinct first indices 
    Is:= {$1..nIs}, #set of first indices
    Pairs1:= combinat:-choose(Is,2),

    #fundamental set of lists (indexed variables without a stem):
    i, j, parms:= {seq(seq([i,j], j= 0..9), i= Is)}, 

    permutations_by_index:= [
        table([2= 3, 3= 2]), #permutations of 1st index
        table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]) #... of 2nd index
    ],

    #Exclusion conditions:
    #---------------------
    #Build a table of sets of the 2nd index grouped by their 1st index:
    ClassifyI:= proc(S)
    option cache;
        (op~)~(2, ListTools:-Classify(x-> x[1], S))
    end proc,

    op2:= proc(S) option cache; op~(2,S) end proc,
    
    TS_5t:= table([0=(), 1=1, 2=2, 3=3, 4=1, 5=(1,2), 6=(1,3), 7=2, 8=(2,3), 9=3]),
     TS_5:= proc(j) option remember; TS_5t[j] end proc,
     TS_6:= table([{1,2}={0,1,2}, {1,3}={0,1,3}, {2,3}={0,2,3}]),

    Condition:= table([
        1= proc(S) option cache; op~(1,S) = Is end proc,
        2= (S-> max(op2(S)) >= 4),
        3= (S-> nops({entries(ClassifyI(S))}) = nIs),
        4= (S->
            Condition[1](S) and
            not ormap(k-> ClassifyI(S)[k] subset [{0,1,4}, {0,2,7}, {0,3,9}][k], Is)
           ),
           5= (S-> TS_5~({op2(S)[]}) = Is),
           6= (S->
               Condition[1](S) and
               not ormap(p-> `union`(map2(index, ClassifyI(S), p)[]) subset TS_6[p], Pairs1)
              )
       ]),
    Conditions:= [1,2,5,4,6], #The order is for efficiency.
            
    #Decide whether a permutation's orbit is allowed to be represented:
    AllowOrbit:= S-> andmap(k-> Condition[k](S), Conditions),
    
    #=================
    #End of setup code

    #Set permutations to identity for unmentioned indices:
    T:= proc(T::table, k) option remember; `if`(assigned(T[k]), T[k], k) end proc,

    #Combine individual indices' permutations into a permutation function for lists:
    conds:= S-> map(x-> T~(permutations_by_index, x), S),

    #For any tuple of lists, compute its orbit under `conds`. Return a
    #representative of the orbit.
    Orbit:= proc(S)
    local r:= S, R, ao;
        do 
            if AllowOrbit(r) then R[r]:= () else return fi 
        until assigned(R[(r:= conds(r))]);        
        {indices(R, 'nolist')}[1] #Return lexicographic min as representative.
    end proc,

    Combos, #:=It:-Combination(nops(parms), tuple_size)

    DoChunk:= (rn::[posint,posint])->  #starting rank and number of iterates
        local 
            Combo_:= Object(Combos, 'rank'= rn[1]),
            has:= ModuleIterator(Combo_)[1], C:= output(Combo_)
        ;       
        (to rn[2] while has() do Orbit(parms[[seq(C+~1)]]) od),
        
    ModuleApply:= proc(tuple_size::posint, {sequential::truefalse:= false})
    local it;
        Combos:= It:-Combination(nops(parms), tuple_size);
        for it to 2 do #kludge to workaround error on 1st run after restart
            try
                return
                    `if`(sequential, map, Threads:-Map['tasksize'= 1])(
                        DoChunk,            
                        {It:-SplitRanks(Number(Combos))[]}
                    )
            catch:
                if it=1 then
                    printf(
                        "Error trapped, "
                        "so efficiency measures are invalid for this run.\n"
                    )
                else
                    error
                fi
            end try
        od
    end proc 
;
end module:

Please ask for definitions of any terms that you're not familiar with, especially those that are italicized.

From a mathematical perspective, tables and functions operating on finite sets are the same thing. This is probably obvious to you. It's only other nonmathematical factors that determine whether they're coded as tables or procedures. And a procedure with option remember or option cache is a hybrid of a table and a procedure.

Note that the tables in permutations_by_index, after being adjusted in T, are not merely functions of {$1..3} and {$0..9}; they are permutations (i.e., bijections with the same domain and range). I assume that this will always be the case---that you'd never have a similiar problem where these tables were not permutations. This observation is key to the rest of this article. Also assume for the rest of this article that the permutations operate on finite sets.

Theorem 1: The Cartesian product of permutations is also a permutaion, operating on the Cartesian product of the sets that the factor permutations operate on.
Proof: It's obvious; or you could likely find a proof in most combinatorics textbooks.

In our code, conds is precisely the Cartesian product of permutations_by_index, as adjusted by T.

Theorem 2: Let P be permutation on A. Let kA be the set of all k-combinations of A. Then the extension of to kA is also a permutation. In particular, applied to any k-combination is a k-combination.
Proof: I leave it for the reader.

Definition 1: Let P be a permutation and S a member of its domain. The orbit of S is the set {S, P(S), P(P(S)), ...}, which is necessarily finite. (This is a standard definition.)

Theorem 3: The orbits of a permutation defined on form a partition of A.
Proof: It's fairly obvious, and you could likely find it in numerous books.

Understanding all of the above is fundamental to understanding our program.

Definition 2: Let be a boolean predicate defined on a set A and let P be a permutation of A. We'll call C homorbital if it being true for one member of an orbit implies that it's true for all members. We call it necessarily homorbital if it's homorbital for any permutation. (These are my own definitions.) 

Now we're ready to address our specific conditions. Let A and B be sets with permutations PA and PB, and let P be the natural permutation defined on the set SS of k-combinations of the Cartesian product AxB. Let C1 be the condition defined on SS by C1(S) is true iff the set of first elements of S is all of A. (This is our condition 1). Then C1 is necessarily homorbital.
Proof: PA(A) = A for any permutation PA, by the definition of permutation.

With the same setup, let C2 be the condition that the maximum of the second elements of S is greater than or equal to 4. Then C2 is homorbital for our specific permutations.
Proof: Our permutation of the second position doesn't exchange any element greater than or equal to 4 with one less than 4.

Our condition 3 is necessarily homorbital. I won't bother proving this because you're no longer using it.

Conditions 4 and 5 are not homorbital, but you said that their being false necessarily disqualifies their counterparts, which amounts to the same thing as being homorbital for the purpose of our program.

What about condition 6? Does it also disqualify the counterpart(s)?

 

@emendes Variable C is defined in condition 3 and used in condition 4. Is condition 3 no longer needed?

Here is the "Exclusion conditions" section of code where I've added your condition 5 and commented out condition 3. The order of the conditions is the fastest order to evaluate them in.

	#Exclusion conditions:
	#---------------------
	#Build a table of sets grouped by their 1st index:
	ClassifyI:= proc(S)
	option cache;
		ListTools:-Classify(x-> x[1], S)
	end proc,
	
	TSt:= table([0=(), 1=1, 2=2, 3=3, 4=1, 5=(1,2), 6=(1,3), 7=2, 8=(2,3), 9=3]),
        TS:= proc(j) option remember; TSt[j] end proc,

	#Decide whether a permutation's orbit is allowed to be represented:
	AllowOrbit:= S->
		local J;
		op~(1,S) = Is					    #condition 1
			and
		max((J:= op~(2,S))) >= 4                      	    #condition 2
			and
		TS~(J) = Is                                  	    #condition 5
		(* Comment out condition 3.
			and	
		nops({entries((op~)~(2,(C:= ClassifyI(S))))}) = nIs #condition 3
		*)
			and not                                     #condition 4:
		ormap(k-> op~(2,ClassifyI(S)[k]) subset [{0,1,4}, {0,2,7}, {0,3,9}][k], Is),

Change notes:

  1. ClassifyI should have option cache if condition 3 is removed.
  2. Anything (including line breaks) between (* and *) is a comment. I use this style for commenting out code.
  3. Change C to ClassifyI(S) in condition 4.

Do you want condition 5 to have the same property that you wanted for condition 4, namely that it being false automatically disqualifies the counterpart of S? That's the way that's coded above.

 

@Joe Riel @emendes

Hi Joe,

We need some help with Iterator. We have a very weird bug that's OS dependent, and in the case of Windows or MacIntosh, it only occurs on the first run after a restart. In the case of Linux, it occurs every run. I suspect that it's related to compilation of a rank-split Iterator used under Threads.

On the Linux machine only, the error occurs while trying to Unrank 1470301, which is the first value in the 27th of 36 pairs produced by SplitRanks. On Mac and Windows 10, the error is about the nonexistence of some weirdly (randomly?) named file.

Best regards,
Carl


Eduardo,

Please add option compile= false to the It:-Combination command and test again on all available OS's.

Second thing to test, independently, and only on the server: Add the option numtasks= 26 to the It:-SplitRanks command.

Third thing to test: Use option sequential in the call to OrbitPartition.

Another significant difference between sets and lists is that membership checks in large sets are faster than those in long lists because it makes use of the fact that the set entries are stored sorted. In addition to the differences mentioned by Acer, I think that this is the only other significant difference.

@emendes I had similar weird errors the 1st (and only the 1st after each restart) time that I tried to run it on each of my two computers: one error coming from Iterator about a random-seeming particular iterate, and one about some randomly-generated file name not existing. Funny thing, though: Both of my computers are Windows 10!! So, please just try again. I have a vague suspicion that the issue is related to Iterator compiling its results (perhaps in separate threads).

The box that contains the module code is called a code edit region. You can use the Insert menu to insert one anywhere in a worksheet. To "load" the code in the region, right click on it and select "Execute Code" from the context menu that pops up. If the code has errors, you'll be notified, but you may need a magnifying glass to locate them precisely: Their approximate locations in the code are marked by red squiggles the same size as the red fibers in US currency.

@Kitonum Here is some cleaner and faster code that does the same thing. My only goal was to make it cleaner; the faster came as an unexpected bonus.

restart:
A:= <1,2,2; 2,1,2; 2,2,3>:
T:= <L[1], -L[2], L[3]; L[1], L[2], L[3]; -L[1], L[2], L[3]>:
NextTriples:= subs(R= convert(T.A, listlist)[], L-> R):

T1:= CodeTools:-Usage(seq((NextTriples~@@n)([[3,4,5]])[], n= 1..10)):
memory used=24.48MiB, alloc change=46.01MiB, cpu time=187.00ms,
real time=189.00ms, gc time=46.88ms

nops([T1]);
                             88572

#Compare with:
NextTriple:=proc(L)
local i;
[seq(op([[L[i,1]-2*L[i,2]+2*L[i,3], 2*L[i,1]-L[i,2]+2*L[i,3], 2*L[i,1]-2*L[i,2]+3*L[i,3]], [L[i,1]+2*L[i,2]+2*L[i,3], 2*L[i,1]+L[i,2]+2*L[i,3], 2*L[i,1]+2*L[i,2]+3*L[i,3]],
[-L[i,1]+2*L[i,2]+2*L[i,3], -2*L[i,1]+L[i,2]+2*L[i,3], -2*L[i,1]+2*L[i,2]+3*L[i,3]]]), i=1..nops(L))];
end proc:
T2:= CodeTools:-Usage(seq((NextTriple@@n)([[3,4,5]])[], n= 1..10)):
memory used=37.53MiB, alloc change=-6.01MiB, cpu time=485.00ms, 
real time=495.00ms, gc time=140.62ms

evalb([T1]=[T2]);
                              true

 

@emendes Garbage collection (gc) happens automatically. Every time that you see either number (memory or time) change in the status bar (bottom right of your worksheet window), there has been a gc. With the code that we're working on, you'll notice that that happens about once per second. There's no danger in doing your own gc, but it's unlikely that you'll improve upon the efficiency of the automatic ones, and there's a good chance that you'll decrease the efficiency. That's why it's discouraged.

@JAMET 

There are many ways to do it. The way that I would recommend would depend largely on two conditions:

  1. Do you want to use a very large n, say, n > 1000?
  2. Do you want to return the entire sequence of iterates, or just its final entry?

The commands foldl and foldr can be used to compose a binary operator with itself. If the operator is associative (as the present case is), these commands give the same result. For example,

foldr(`&*`, [3,4,5] $ 9)

The procedure that you were trying to write could be written like this:

Tri:= proc(Y::And(list(algebraic), 3 &under nops), n::posint)
local r:= Y;
    to n-1 do r:= Y &* r od;
    r
end proc:

You shouldn't declare a variable global just because you use it in a procedure. That's only needed if you make an assignment to it.

@emendes If your attempt to use Grid is pursuant to the same problem that we're working on in the other thread, then don't bother. I just finished parallelizing our OrbitPartition with Threads, which is much more memory efficient than Grid (if you can satisfy its rather strict requirements).

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