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


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)))));

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:


gc(); #just to refresh status bar


OrbitPartition:= module()





(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}]),
            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



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





@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. 



Digits := trunc(evalhf(Digits));



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)))]:


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))})


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]


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)[]





@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.


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.

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):

        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);


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])





@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.

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