Question: How to tweak a procedure to avoid choking as much as possible when using Grid?

Hi,

I need to run the following procedure a couple of million times. Although it works, Maple sometimes chokes for no apparent reason (if there is a reason, please let me know). I was wondering whether an expert could help me tweak the procedure (or possibly rewrite it) to achieve the best possible performance. I am planning to use Grid:-Map or, if possible, Threads:-Map.

generateNonlinearModelsPlus := proc(model::list,fullmodel::list,vars::list:=[x,y,z])
description "This function generates a list of all models with one more monomial from the full model":
local tab::table(),n:=nops(model),i,j,k:=1,ans,terms,aaa,allmoncoefThreads:
# local procedure
allmoncoefThreads := proc(f::list,vars::list)
description "This function finds the monomials multipled by their coefficients for each expression (equation) of a list.":
local n:=numelems(f),i,mon:=[seq](0,i=1..n),M,cc:=[seq](0,i=1..n),ans:
for i from 1 to n do
  cc[i]:=[coeffs](expand(f[i]),vars, 'M'):
  mon[i]:=[M]:
end do:
ans:=[seq](zip((ww,vv)->ww*vv,cc[i],mon[i]),i=1..n):
return(ans)
end proc:
# main part
ans:=zip((w,v)->expand(simplify(v-w)),model,fullmodel): # Find the monomials that are not in model
terms:=allmoncoefThreads(ans,vars): # Separate the monomials
#
for i from 1 to n do
   aaa:=model:
   for j from 1 to nops(terms[i]) do
       aaa[i]:=model[i]+terms[i,j]:
       tab[k]:=aaa:
       k:=k+1:
   end do:
end do:
tab:=convert(tab,list):
return(tab):
end proc:

Here is an example of how I run it: 

model:=[y*alpha[1, 2], z*alpha[2, 3], x^3*alpha[3, 10] + x*alpha[3, 1] + alpha[3, 0]]:

fullmodel:=[x^3*alpha[1, 10] + x^2*y*alpha[1, 11] + x^2*z*alpha[1, 12] + x*y^2*alpha[1, 13] + x*y*z*alpha[1, 14] + x*z^2*alpha[1, 15] + y^3*alpha[1, 16] + y^2*z*alpha[1, 17] + y*z^2*alpha[1, 18] + z^3*alpha[1, 19] + x^2*alpha[1, 4] + x*y*alpha[1, 5] + x*z*alpha[1, 6] + y^2*alpha[1, 7] + y*z*alpha[1, 8] + z^2*alpha[1, 9] + x*alpha[1, 1] + y*alpha[1, 2] + z*alpha[1, 3] + alpha[1, 0], x^3*alpha[2, 10] + x^2*y*alpha[2, 11] + x^2*z*alpha[2, 12] + x*y^2*alpha[2, 13] + x*y*z*alpha[2, 14] + x*z^2*alpha[2, 15] + y^3*alpha[2, 16] + y^2*z*alpha[2, 17] + y*z^2*alpha[2, 18] + z^3*alpha[2, 19] + x^2*alpha[2, 4] + x*y*alpha[2, 5] + x*z*alpha[2, 6] + y^2*alpha[2, 7] + y*z*alpha[2, 8] + z^2*alpha[2, 9] + x*alpha[2, 1] + y*alpha[2, 2] + z*alpha[2, 3] + alpha[2, 0], x^3*alpha[3, 10] + x^2*y*alpha[3, 11] + x^2*z*alpha[3, 12] + x*y^2*alpha[3, 13] + x*y*z*alpha[3, 14] + x*z^2*alpha[3, 15] + y^3*alpha[3, 16] + y^2*z*alpha[3, 17] + y*z^2*alpha[3, 18] + z^3*alpha[3, 19] + x^2*alpha[3, 4] + x*y*alpha[3, 5] + x*z*alpha[3, 6] + y^2*alpha[3, 7] + y*z*alpha[3, 8] + z^2*alpha[3, 9] + x*alpha[3, 1] + y*alpha[3, 2] + z*alpha[3, 3] + alpha[3, 0]]:

vars:=[x,y,z]:

ans:=generateNonlinearModelsPlus(model,fullmodel,vars)

Many thanks.