Carl Love

Carl Love

28025 Reputation

25 Badges

12 years, 313 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@acer I'm sure that I was the original author of that code, but it looks like it's been put through the wringer a few times, extracting all comments and formatting, and incorrectly modified a bit, and converted to 2D Input and then back to plaintext.

@emendes 

Re: Homorbital: As I said, that one is my own definition. Perhaps the concept is defined elsewhere with a different name,

Re: counterparts: Don't think in terms of what is removed. Nothing is removed in this program that uses an Iterator. Rather, only items that are desired for the final output are included in the first place. That's why this program is memory conservative. The more conditions that you add, the more time it will use, but the less memory.

Re: arguments: Yes, I'll add those things as arguments.

You say that parms can change. How exactly? Will there always be two indices? Is the first index always $1..nIs?

Can you think of names for the 6 conditions? It would be useful to have key names to replace the 1-6 used as table indices.

Can you describe how the tables of sets for conditions 4-6 are created? It would be useful to create them programmatically from some higher concept if possible.

Here is a worksheet with all-new and faster code, including passing the arguments. I do not pass nIs because it can be derived from parms.

Symmetries.mw

Please test scenarios that gave errors before. I might have fixed some of them. I can't test because I don't have Linux or a machine with 36 processors.

@nm Maple does have a limited ability to translate some Mathematica. See ?MmaTranslator.

@crashoverwide Oops, I have any intermittent GUI bug that causes minus signs to show up as K, even on graph axes. I just automatically translate them in my mind; unfortunately that's how common this bug is. So I just realized that I posted such a graph above. So, here's the graph with regular minus signs.

Does anyone know anything about this bug. It happens on both my computers. Plus signs show as C, and there are many other weird transliterations.

@dingtianlidi Here's a bonus problem, since you seem interested. Using the same type of graphical analysis, we can convert the same integration region to polar coordinates. For this, consider a ray emanating from the origin. Through which curve does it enter the region? Through which curve does it exit the region?
 

f1:= sqrt(4*x-x^2): f2:= 2*sqrt(x):

plots:-display(
    plots:-shadebetween(
        piecewise(x<2, f1, x), f2, x= 0..4,
        color= magenta
    ),
    plots:-shadebetween(f1, x, x= 2..4, color= cyan),
    plot(x, x= 0..4, color= black, linestyle= dash, thickness= 2),
    gridlines= false
);

Polar:= f-> solve(eval(f, [x,y]=~ r*~[cos,sin](theta)), r):

g1:= Polar(y=f1);

0, 4*cos(theta)

g2:= Polar(y=f2);

0, 4*cos(theta)/sin(theta)^2

g3:= Polar(x=4);

4/cos(theta)

Int(r, [r,theta]=~ [g1[2]..g3, arctan(0,4)..arctan(4,4)]) +   #cyan
Int(r, [r,theta]=~ [g1[2]..g2[2], arctan(4,4)..arctan(4,0)]); #magenta

Int(r, [r = 4*cos(theta) .. 4/cos(theta), theta = 0 .. (1/4)*Pi])+Int(r, [r = 4*cos(theta) .. 4*cos(theta)/sin(theta)^2, theta = (1/4)*Pi .. (1/2)*Pi])

value(%);

32/3-2*Pi


 

Download Polar_conversion.mw

@Steve Roper It's safer to generalize that if-statement to 

if not [args]::list(numeric) then return 'procname'(args) fi

Thus it'll work correctly regardless of the procedure's parameter names.

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.

First 198 199 200 201 202 203 204 Last Page 200 of 708