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

@Carl Love For your work today, please use this newest and fastest version: Symmetries.mw

@Tyttemus It's mathematically impossible to do what you want.

Looking through the help pages (in Maple 2020) whose names begin "updates," the oldest that I found was for Maple 4.0. At least that's something older than Maple 6.

@acer I failed to consider that in my definitions, although I was aware of that in the back of my mind! Thanks for pointing that out, and I'm glad to know that you're still reading this thread and monitoring the accuracy of what I say.

Fortunately, for the situation in this program, we are dealing only with sets of sets of lists of integers (type set(set(list(integer)))), and certainly all that I said applies to those. Agree? Would you be able to modify my definitions to encompass your example? Or is my whole thesis irreparably flawed?

@emendes 

1. Feasibility: Since evalf(7*binomial(60,7)) = 2.7*10^9, I'd guess that that's feasible, but near the upper limit of feasibility.

2. nterms: Yes, it totally makes sense. It's standard combinatorics.

7. Keep me updated on any need for new conditions. And be wary that I'm not sure that I've implemented the conditions correctly.

8. One problem at a time: Suppose that you wanted to use CondCheck with two (or more) different configurations of the index permutations and\or condition tables. Using a "basic" module (which is what we have), you need to re-call OrbitPartition every time that you want to change the configuration. With an object module, you can set different configurations and assign them to names, and those "objects" will exist simultaneously in memory, so that you can switch configurations for CondCheck just by using the different names.

9. I'd like to know the "parallel efficiency", which I define as  "cpu time" / ("real time" * "number of processors").

 

@emendes It's an excellent question, really the key to the whole working of the program, and fundamental to the working of Maple itself as a whole. Consider the return-value line from Orbit:

    {indices(R, 'nolist')}[1] #Return lexicographic min as representative.

At the time that this line is executed, table contains (as its indices) all of the counterparts of as well as S itself (i.e., the orbit of S). (In our case, there can only be one counterpart, but this code will work for any number of counterparts, i.e., any orbit length.) The {...} forms those indices into a set. (The nolist avoids adding extra [...], which are not needed here and just contribute to garbage.) Now notice the [1]: That gets the first element of the set as the return value. Now, you may ask, "How do you define the first element of a set?" Given any two immutable[*1] objects and y in Maple, it has an ordering scheme that is consistent across sessions[*2], even across different computers, such that if x precedes y, then x will come before y in any set whatsoever. It's an absolutely brilliant piece of design on the part of the Maple kernel authors. Thus, I know that that [1] will always select the desired member of the orbit, and thus leave out what you called S'. 

[*1] Immutable objects: For the purpose of this post, immutable objects are those whose "identity" is based solely (by design!) on "content" rather than by memory address. Immutable objects include all numbers, names, algebraic expressions, and lists, sets, and sequences whose elements are themselves immutable. Objects that are not immutable (aka mutable) include tables and all types of arrays (Vectors, Matrices, Arrays, rtables).

[*2] For the purpose of rigorous testing of programs (and pretty much only for that purpose), it's possible to choose different ordering schemes. I think that there are 8 choices. See help page ?maple.

 

 

@emendes 

1. Problem size: I forgot to include the tuple size as a factor in the proportionality estimates of the problem size (both run time and memory consumption). The p is proportional to t*binomial(n*binomial(n+d,d), t), where t is tuple size, n is the number of 1st indices, and d is the degree. For the current (30,15) problem, we have evalf(15*binomial(30,15)) ~ 2.3*10^9. For the (60,30) problem, it's evalf(30*binomial(60,30)) ~ 3.5*10^18. So the (60,30) problem is more than a billion times larger than the (30,15) problem. Dividing the problem into products of 3 binomial(20, ...)s only helps a tiny amount because

evalf(30*add((mul@binomial~)~(20, combinat:-conjpart~(combinat:-partition(30,3)))));
             8.7*10^17

Thus, I think that the (60,30) problem is totally infeasible.

2. nterms and FullDeg: Okay, I totally understand it now, and I've modified the code accordingly.

3. Combining conditions: It's condition 4, not condition 5, that can be combined with condition 6; and I've modified the code accordingly.


Code changes:

1. Conditon names: 
        condition 1  ==>  AllI (as before)
        condition 2  ==>  FullDeg
        
condition 3  ==>  removed
        condition 4  ==>  incorporated into FullDim (condition 6)
        condition 5  ==>  ValidJ
        
condition 6  ==>  FullDim

2. The cut-off value for FullDeg is named FullDegJ, and it's computed from nterms in Init as you specified.

3. The singleton subsets of Is are now included in SubsetsI.

4. The module now exports a procedure named CondCheck that does what you asked for in an email.

5. Since the problem must be first "set up" before CondCheck can be run, there's now a truefalse module local named Setup.

6. Since Setup must start out false before any problem is run there's a procedure ModuleInit that makes it so.

7. Since you might want to set up a problem without actually iterating through tuples just so that you can use CondCheck, I included a mechanism for that: Just specify a tuple size of 1. You'll still be able to call CondCheck with a tuple of any size.

8. The module can only be "set up" for one problem at a time. If you need to change that, let me know. It would require changing the module to an object module.

9. I improved the runtime by a factor of 3 - 4 (based on very limited testing) by an extremely subtle manipulation of cache/remember tables inside of cache/remember tables. Discussing the precise details of how this works would only be suitable for the most-expert-level course in Maple programming, so I'll leave it at that. You can see the details in the very short procedures JsbyI and CT.

10. Cache sizes are set based on the number of processors.


I'll post the new code in a new Answer because this subthread is getting long.

 

@emendes Thanks, I'm beginning to understand the condition Nonlinear. Perhaps it should be renamed FullDeg?

What is nterms? Have you coded it? Or should it be an input parameter?

I suspect that conditions 5 and 6 can be combined, as you suggest. I'll need to look at it further. If so, I think there'll be an improvement in run time. Checking the conditions is the bottleneck for for run time.

binomial(30,15) = 155,117,520. Divided over 36 processors, that's 4,308,820. The real run time should be proportional to that, and 4.3M seems easily doable.

BUT, memory needed will be proportional the 155M. Keep an eye on that memory number in the status bar, or use an OS tool to monitor memory usage. We don't want it to start swapping; it's very slow! If swapping is needed, we can manage it from within the program by rolling off large chunks of the processed elements to disk.

@emendes From when you first introduced condition 2 (now called Nonlinear), the critical value was 4. Now you're saying that it's 10?? That's impossible, because the maximum value of the 2nd index is 9 (j= 0..9).

We have:

Terms:            1, x, y, z, nonlinears...
Corresponding j:  0, 1, 2, 3, 4, ...

All terms after 3 are nonlinear (in the 3D case). I want to generalize that to "all terms after nIs are nonlinear".

In other news:

  1. I changed procedure name conds to Symm.
  2. I changed condition name NoPairs to FulDeg (short for full degree),
  3. I changed set name Pairs1 to SubsetsI, and it contains all subsets of sizes 2..nIs-1.

Is this all good?

@emendes I am asking whether 4 in Nonlinear can be replaced by nIs+1. I'm not asking what max(Js(x)) is equal to.

Linux: Run it again. See if it takes 90 minutes again. I assume that you mean "real time" > 90?

@emendes I meant that you could suggest areas for commenting, and I would write them, assuming that they were areas that you have questions about.

Am I correct in assuming that the 4 in the NonLinear condition is actually nIs + 1?

I will extend the NoPair concept as needed.

The PDF files are password locked. If you wish to keep it so, you could email me the password,

@emendes Okay, I've made several changes. Highlights:

  1. All conditions have been recoded to be indexed by names. The names that I (provisionally) used are AllINonlinearDifferent`#I-D:I``#I-D:J`NoPairs, and Symm. I'm not totally pleased with the 4th and 5th of these. Can you suggest something else, perhaps names that won't need quotes? It's fine if you use technical language.
  2. Thank you for code snippet with nIs and degree. I made them input parameters, and I removed parms, which is now constructed internally.
  3. I give some more control over the bugs associated with compiled and\or multithreaded Iterators, I've added keyword parameters numtasks and nocompile to the sequential that I already had. The default value of numtasks is kernelopts(numcpus), which is I believe 36 on your Linux machine. So, please try numtasks=35 for the example that errored out. If it still errors out, try adding nocompile also. That'll be a little slower, but it'll give @Joe Riel more info for when he gets around to looking at this. (I had no trouble doing the 8-tuples on my puny machine, but it takes 90 minutes.)
  4. Several variables have been given more logical names. Feel free to suggest other such changes. I do strongly prefer a short name with a comment over a longer name that is itself the comment.
  5. I've added several more comments, and elaborated several existing comments. If there are areas for which you'd like more comments, please suggest.

Please attach that PDF to your Reply.

I am eager to code the 7th condition: symmetry.

I think that you should mark this Answer (the parent of this Reply) as the Best Answer.

Here is the new work:

restart:

gc(); #just to refresh status bar

 

OrbitPartition:= module()

 

 

kernelopts(numcpus);

4

(nIs, deg):= (3,2) #number of 2nd indices will be binomial(3+2,2) = 10.
:
permutations_by_index:= [
    table([2= 3, 3= 2]),
    table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7])
]:
CondTables:= Record(
    `#I-D:I`= table([1= {0,1,4}, 2= {0,2,7}, 3= {0,3,9}]),
    `#I-D:J`=
        table([
            0=(), 1=1, 2=2, 3=3, 4=1, 5=(1,2),
            6=(1,3), 7=2, 8=(2,3), 9=3
        ]),
    "NoPairs"= table([{1,2}={0,1,2}, {1,3}={0,1,3}, {2,3}={0,2,3}])
):
Args:= nIs, deg, permutations_by_index, CondTables
:

Result1:= CodeTools:-Usage(OrbitPartition(8, Args)):

memory used=0.60TiB, alloc change=0.92GiB, cpu time=3.37h, real time=95.78m, gc time=7.41h

nops(Result1);

2161507

Result2:= CodeTools:-Usage(OrbitPartition(4, Args, numtasks= 5)):

memory used=0.79GiB, alloc change=0 bytes, cpu time=6.05s, real time=2.31s, gc time=1.64s

nops(Result2);

2296

 

Download Symmetries.mw

@Grigoriy Yashin Yes, trig is one of the simplification algorithms that combine knows. To apply all of the algorithms that it knows, just remove the word trig. See ?combine, which lists the names of the algorithms. You may also want to apply simplify after combine.

@vv Thanks. Now I understand option threadsafe. It's just a way to flag a procedure as threadsafe in case the kernel incorrectly suspects otherwise; it overrides that suspicion.

Google has a much superior search algorithm, which is applicable here because Maple help pages are Google indexed (since they're somewhat duplicated as regular web pages). 

It burns me up that MaplePrimes Answers are not properly searchable by Google! The best hit that you'll get is the Answers page of the user who wrote the Answer. Since I have about 4000 Answers, someone Google searching for one of my Answers only knows that they have a positive hit; they'd need to wade through all the Answers to find it.

I got your email about this problem, where you attached a worksheet, which I re-attach here.

No Maple code that I have ever written (or ever will write) is intended to be used by Maple's 2D Input. Considering the length of the code, I refuse to look for the error in the 2D Input. If you redo it in the 1D Input form (code in a reddish monospaced font) that I originally posted it in, then I'll look at it. 


 

restart

Digits := trunc(evalhf(Digits));

15

(1)

ODEs := [diff(f(eta), `$`(eta, 3))+A^2+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2+beta*((diff(g(eta), eta))^2-g(eta)*(diff(g(eta), `$`(eta, 2)))-1), lambda*(diff(g(eta), `$`(eta, 3)))+(diff(g(eta), `$`(eta, 2)))*f(eta)-g(eta)*(diff(f(eta), `$`(eta, 2)))]:

`<,>`(ODEs[])

Vector(2, {(1) = diff(f(eta), eta, eta, eta)+A^2+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2+beta*((diff(g(eta), eta))^2-g(eta)*(diff(g(eta), eta, eta))-1), (2) = lambda*(diff(g(eta), eta, eta, eta))+(diff(g(eta), eta, eta))*f(eta)-g(eta)*(diff(f(eta), eta, eta))})

(2)

LB, UB := 0, 1:

BCs := [`~`[`=`](([D(f), f, g, (D@@2)(g)])(LB), [1+B1*((D@@2)(f))(0), 0, 0, 0])[], `~`[`=`](([D(f), D(g)])(UB), [A, 0])[]];

[(D(f))(0) = 1+B1*((D@@2)(f))(0), f(0) = 0, g(0) = 0, ((D@@2)(g))(0) = 0, (D(f))(1) = A, (D(g))(1) = 0]

(3)

Params := Record(A = .9, B1 = .5, beta = .5, lambda = .5):

NBVs := [-((D@@2)(f))(1) = `C*__f`]:

Cf := `C*__f`:

Solve := module () local nbvs_rhs, Sol, ModuleApply, AccumData, ModuleLoad; export SavedData, Pos, Init;  nbvs_rhs := `~`[rhs](:-NBVs); ModuleApply := subs(_Sys = {:-BCs[], :-NBVs[], :-ODEs[]}, proc ({ A::realcons := Params:-A, B1::realcons := Params:-B1, beta::realcons := Params:-beta, lambda::realcons := Params:-lambda }) Sol := dsolve(_Sys, _rest, numeric); AccumData(Sol, {_options}); Sol end proc); AccumData := proc (Sol::{Matrix, procedure, list({name, function} = procedure)}, params::(set(name = realcons))) local n, nbvs; if Sol::Matrix then nbvs := seq(n = Sol[2, 1][1, Pos(n)], n = nbvs_rhs) else nbvs := `~`[`=`](nbvs_rhs, eval(nbvs_rhs, Sol(:-LB)))[] end if; SavedData[params] := Record[packed](params[], nbvs) end proc; ModuleLoad := eval(Init); Init := proc () Pos := proc (n::name) local p; option remember; member(n, Sol[1, 1], 'p'); p end proc; SavedData := table(); return  end proc; ModuleLoad() end module:

colseq := [red, green, blue, brown]:

Pc := B1 = .5, A = .1, beta = .5:

Ps := [B1 = .5, A = .1, beta = .5]:

Pv := [lambda = [.5, 1, 1.5, 2]]:

for i to nops(Ps) do plots:-display([seq(plots:-odeplot(Solve(lhs(Pv[i]) = rhs(Pv[i])[j], Ps[i][], Pc), [eta, theta(eta)], 'color' = colseq[j], 'legend' = [lhs(Pv[i]) = rhs(Pv[i])[j]]), j = 1 .. nops(rhs(Pv[i])))], 'axes' = 'boxed', 'gridlines' = false, 'labelfont' = ['TIMES', 'BOLDOBLIQUE', 16], 'caption' = nprintf(cat(`$`("\n%a = %4.2f, ", nops(Ps[i])-1), "%a = %4.2f\n\n"), `~`[lhs, rhs](Ps[i])[]), 'captionfont' = ['TIMES', 16]) end do

Error, (in dsolve/numeric/process_input) invalid argument: (B1 = .5)[]

 

``


 

Download 23.mw

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