@acer 

I did try to increase it using kernelopts(stacklimit=...), but on macOS there seems to be a hard barrier around 65xxx, so I couldn't push it any further.

I used Grid:-Map because I wasn't sure whether indets was thread-safe in my original implementation. Since Grid:-Map launches separate Maple processes, I knew that wouldn't be an issue there.

Now I'd like to go back and try the fastest implementation, using Threads:-Map with indets, and see whether indets turns out to be safe in this particular application. If it is, that would clearly be the preferred approach.

@acer 

Thank you very much for these tests and for the alternative implementations. The timing results are extremely helpful.

To answer your question in the postscript: no, the same coefficient alpha[i,j] cannot appear twice within a given model. For each fixed pair [i,j], alpha[i,j] is unique within that model. Therefore, an expression such as

x*z^2*alpha[3,19] + z^3*alpha[3,19]

cannot occur in my actual data.

When running ex53_7.mw on my Mac I got the following results: (I don't know what you did to get that nice look when listing the results).
 

# mine
time[real]( map(w->map(P,eval({w[]},H)),models) );
                            128.024

# mine, Thr
time[real]( Threads:-Map(w->map(P,eval({w[]},H)),models) );
                             89.740

# OP original
time[real]( map(w->map(v->[op](v),indets(w,name) minus svars),models) );
                             55.953

# my tweak of OP original
time[real]( map(w->map([op],indets(w,indexed)),models) );
                             42.986

time[real]( Threads:-Map(w->map([op],indets(w,indexed)),models) );
                             39.301

Please notice that there is no much difference between using Threads or not.  

Finally, I tried running the fastest implementation on the eight-term models, but it ended with a "stack limit exceeded" error on my Mac.

Before posting my question on MaplePrimes, I had already run the Grid:-Map implementation on a Linux machine with 128 CPUs and 2 TB of RAM, and it completed successfully after approximately six hours.

Once I have access to that machine again, I will try the fastest implementation there and report back the results.

@acer 

I tried your third alternative implementation using the complete four-term dataset (45,770 models). For some models, however, it returns only three elements instead of the expected four. 

It fails for monomials x^a*y^b*z^c*alpha[i,j] with a=b=c=0.  

@acer 

Thank you. I now understand why the precise structure of the data is important.

There is a very large number of modestly sized lists in models; each model is always a list of exactly three polynomial expressions, corresponding to a three-dimensional autonomous system.

All models have exactly the same number of monomial terms. In the example I posted, each model has four terms in total. In the actual computation, every model has eight terms in total, and there are more than 169 million such models.

The only non-indexed names are

svars := {x,y,z};

The only indexed names are the coefficients alpha[i,j], where i and j are positive integers.

Every monomial has exactly the form

x^a*y^b*z^c*alpha[i,j]

where a, b, and c are nonnegative integers satisfying

a+b+c <= 3.

Thus, every model is a cubic polynomial system. There are no functions, rational expressions, additional symbolic variables, or other indexed names.

My objective is simply to extract, for each model, the set of distinct index pairs [i,j] corresponding to the coefficients alpha[i,j] (Carl's notation). Since the polynomial structure is completely fixed, I was wondering whether there might be a much faster approach than using a general symbolic routine such as indets.

I also agree with your point that optimizing the serial operation before attempting to parallelize it is probably the right strategy.

@acer Many thanks.  Would it be possible for you to provide a small example of the implementation you have in mind, particularly the approach based on tables (or another more suitable data structure)? I would greatly appreciate seeing a concrete implementation, as this would help me better understand how to adapt the procedure for large-scale computations. I do apologize for asking this again and for not having learned from the countless times you and Carl Love have helped me in the past.

@janhardo Thank you.   

@janhardo Pretty impressive. Thank you very much. Unfortunately, I tend to overcomplicate things when writing a procedure for the first time. Also, Vars is no longer needed.

@C_R There is no problem with the uploaded example itself. However, when I run the procedure over a large pool of models (for example, systems with 6 monomials instead of 5, resulting in about 203685 systems in total), Maple eventually stops progressing at some point. At that stage, CPU usage drops from nearly 100% across all cores to almost zero, and Maple appears to stall indefinitely without producing an error message.

I came up with the following simplified procedure that passes the ThreadSafetyCheck,

generateNonlinearModelsPlus_fast := proc(model::list, fullmodel::list, vars::list := [x,y,z])
description "Generate all models with one additional monomial from fullmodel":
local n, i, j, diff, coeff, mon, out, k;

    n := nops(model);
    out := table();
    k := 0;

    for i to n do
        # One expansion is enough; simplify is usually much more expensive here.
        diff := expand(fullmodel[i] - model[i]);

        if diff <> 0 then
            coeff := [coeffs(diff, vars, 'mon')];
            mon := [mon];

            for j to nops(coeff) do
                k := k + 1;
                out[k] := subsop(i = model[i] + coeff[j]*mon[j], model);
            end do;
        end if;
    end do;

    return [seq(out[i], i = 1 .. k)];
end proc:

@sand15 

The idea is to generate all systems containing a total of (n+1) monomials from a pool of systems containing (n) monomials that satisfy a given condition. Each new system is obtained by adding a single monomial selected from the set difference between the full system and the corresponding (n)-monomial system.

As an example, the full system may correspond to a 3D polynomial system containing all possible monomials up to degree 3.

@sand15 Many thanks. Here it is an example of how to run the procedure.

Example_Generate_Models.mw

As for threadsafety, I knew that zip was not threadsafe, thus using Grid.

@acer Thank you for the explanation. fixup saves me a lot of time when transferring my results from Maple to LaTeX.  

@acer Thank you, but I think I may have missed something. fixup returns the same input expression expr, doesn't it? How was it modified?

@acer Sorry to reopen the thread, but do you happen to know how to fix z*x? I mean…

from

[sigma*y, (-sigma - 1)*y + z*x, (rho - 1 - z)*beta - (x + y)*x]

to

[sigma*y, (-sigma - 1)*y + x*z, (rho - 1 - z)*beta - (x + y)*x]

@acer I could not ask for more.  Thank you.  

@acer I could never have come up with such a clean procedure—it's got Carl written all over it. I can't thank him enough for all the help he's given me. And now, you're part of that too.

Example

lorenz := [-sigma*x + sigma*y, rho*x - x*z - y, -beta*z + x*y];
vars := [x,y,z];

map(w->monomialsGB(w,vars),lorenz);

and output is

[[sigma, -sigma], [y, x], [-1, -1, rho], [z*x, y, x], [-beta, 1], [z, y*x]]

The ordering should be [x,y], [x,y,x*z] and [z,x*y].  

Actually I need both versions. I mean Carl's and yours.   

Many thanks.