#3933 storage exhausted in large calculation

[Guillermo P. Ortiz]

Dear Maxima Team,
I am not sure about this issue, but I did not found it in any blog or mailing list.
According to this report a related lisp error with MACSYMA TOP LEVEL memory in CONS
expression is the cause.

++++++++
here the output message:

Maxima encountered a Lisp error:

Condition in MACSYMA-TOP-LEVEL [or a callee]: INTERNAL-SIMPLE-ERROR: The storage for CONS is exhausted. 1634280 pages allocated. Use ALLOCATE to expand the space.

Automatically continuing.
+++++++

and here the bug_report() output:

Please report bugs to:
https://sourceforge.net/p/maxima/bugs
To report a bug, you must have a Sourceforge account.
Please include the following information with your bug report:

---

Maxima version: "5.44.0"
Maxima build date: "2021-04-24 14:52:58"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.12"
User dir: "/home/gortiz/.maxima"
Temp dir: "/tmp"
Object dir: "/home/gortiz/.maxima/binary/5_44_0/gcl/GCL_2_6_12"
Frontend: false

---

Below I copy the code, but to run a couple of tables are required and many pages are displayed. Please, let me know is it is needed anyway. They are interpolated after makelist for eack are doing.
Also I introduce my own sqrt (mysqrt) because I need a different branch cut by default in maxima.

dynamalloc:true;
eq1:ED1*exp(%i*k1*d1)-ED2*exp(%i*k2*d1)+EI1*exp(-%i*k1*d1)-EI2*exp(-%i*k2*d1)=DE0*exp(%i*k*d1);
eq2:k1*ED1*exp(%i*k1*d1)-k2*ED2*exp(%i*k2*d1)-k1*EI1*exp(-%i*k1*d1)+k2*EI2*exp(-%i*k2*d1)=k*DE0*exp(%i*k*d1);
eq3:ED1-ED2*exp(%i*(k2-k)*d)+EI1-EI2*exp(-%i*(k2+k)*d)=DE0;
eq4:k1*ED1-k2*ED2*exp(%i*(k2-k)*d)-k1*EI1+k2*EI2*exp(-%i*(k2+k)*d)=k*DE0;

neq:linsolve([eq1,eq2,eq3,eq4],[ED1,ED2,EI1,EI2]);
E01:4*%pi*q*J0/(%i*c**2*(k1**2-k**2));
E02:4*%pi*q*J0/(%i*c**2*(k2**2-k**2));

neq:subst(DE0=E02-E01,neq);


[ED1,ED2,EI1,EI2]:subst(neq,[ED1,ED2,EI1,EI2]);

sinc(x):=sin(x)/x;
epsM:(k/q)**2*c**2+4*%pi*J0*d/(%i*q*(d1*E01+d2*E02+
    ED1*d1*sinc((k1-k)*d1/2)*exp(%i*(k1-k)*d1/2)+
    EI1*d1*sinc((k1+k)*d1/2)*exp(-%i*(k1+k)*d1/2)+
    ED2*d2*sinc((k2-k)*d2/2)*exp(%i*(k2-k)*d2/2)+
    EI2*d2*sinc((k2+k)*d2/2)*exp(-%i*(k2+k)*d2/2)));

epsM(z):=subst(k=z,epsM);
epsMD2:at(diff(epsM(z),z,2),z=0.001);

mysqrt(x):= if imagpart(sqrt(x))<0 then -sqrt(x) else sqrt(x);

epsM:subst([k1=mysqrt(e1)*q/c,k2=mysqrt(e2)*q/c,d1=80.0,d2=20.0,d=100.0,J0=1.0],epsM(0.001));
s : openr("epsAu_JC.dat");
readline(s);
a:read_nested_list(s);
len:length(a);
are:makelist([a[i][1],a[i][2]],i,1,len);
aim:makelist([a[i][1],a[i][3]],i,1,len);
load("interpol");
epsAu(q):=cspline(are,varname='q)+%i*cspline(aim,varname='q);
close(s);

s1 : openr("epsSiO2.dat");
readline(s1);
readline(s1);
a1:read_nested_list(s1);
len:length(a1);
a1re:makelist([a1[i][1],a1[i][2]],i,1,len);
load("interpol");
epsSiO2(q):=cspline(a1re,varname='q);
close(s1);


epsM:subst([e1=epsSiO2(q),e2=epsAu(q),c=197.33],epsM);
epsM1(x):=subst(q=x,epsM);
epsM2re(x):=float(realpart(epsM1(x)));
epsM2im(x):=float(imagpart(epsM1(x)));

epsMlist:makelist([i,epsM2re(i),epsM2im(i)],i,1.5,3.5,0.04);

thank you in advance for your attention
Regards

´´´

[Robert Dodier]

Ticket moved from /p/maxima/feature-requests/161/

[Robert Dodier]

Guillermo, thanks for the bug report.

My first guess about what's going on is that you are constructing extremely huge expressions in cspline and then trying to do some operations on those which make them even more enormous, and eventually Maxima runs out of memory.

If you attach the data file or files to this bug report, it would help others understand what's going on.

It would also help othes understand if you explain something about what you are trying to do. What is the larger goal you are working towards?

[Guillermo P. Ortiz]

Hello Robert,
I have set dynamalloc to true and a lot enough of ram memory, the maxima process in question only take a half of it, but thinking about your comment I change to a non interpolating version of my code. One issue that I am not sure is how to manage the sqrt redefinition to mysqrt. It is needed because the default definition of the branch cut selected by maxima but not adjustable by user, is not it?
As a mathematical task, briefly I want to solve a 4x4 linear eq. system from there I obtain an expression depending on several parameters and two variables (q,k) and then from here I want to obtain a couple of expression that I made as list and then deliver into two output files.
Regards

[Stavros Macrakis]

Why are you constructing these huge expressions with lots and lots of common subexpressions?
Even the very naive optimize function reduces your epsM from 37k characters to 2650 characters (in string form).

[Guillermo P. Ortiz]

Hello Stavros,
I think that main problem is with mysqrt definition. I have not found any other form to manage the selection of other branch cut for sqrt function defined by Maxima.
I am trying a code version without cspline and seems after 3 hs with core i7 7ma generation processor and 16 GB RAM , that Maxima using 8GB is working to hard to manage it
I will kill it and I could try with optimize before,
thank you in advance
Regards

[Stavros Macrakis]

I can't execute your code beyond the definition of epsM (because of the file load), so I'm not sure what it's doing. But if the goal is numerical (not symbolic) calculation, it is probably unnecessary to carry around those huge symbolic expressions.

A common error in using Maxima is to generate large symbolic expressions which are later evaluated for specific numbers. It is often better just to do the numerical calculation. This is especially true for interpolation.

[Guillermo P. Ortiz]

OK, let me share a code version similar to the above but without cspline. I am totally agree with you about numerical versus symbolic, the problem for me is I am not how to manage the sqrt redefinition according to the neeede branch cut that I said before.

dynamalloc:true;
eq1:ED1*exp(%i*k1*d1)-ED2*exp(%i*k2*d1)+EI1*exp(-%i*k1*d1)-EI2*exp(-%i*k2*d1)=DE0*exp(%i*k*d1);
eq2:k1*ED1*exp(%i*k1*d1)-k2*ED2*exp(%i*k2*d1)-k1*EI1*exp(-%i*k1*d1)+k2*EI2*exp(-%i*k2*d1)=k*DE0*exp(%i*k*d1);
eq3:ED1-ED2*exp(%i*(k2-k)*d)+EI1-EI2*exp(-%i*(k2+k)*d)=DE0;
eq4:k1*ED1-k2*ED2*exp(%i*(k2-k)*d)-k1*EI1+k2*EI2*exp(-%i*(k2+k)*d)=k*DE0;

neq:linsolve([eq1,eq2,eq3,eq4],[ED1,ED2,EI1,EI2]);
E01:4*%pi*q*J0/(%i*c**2*(k1**2-k**2));
E02:4*%pi*q*J0/(%i*c**2*(k2**2-k**2));

neq:subst(DE0=E02-E01,neq);

[ED1,ED2,EI1,EI2]:subst(neq,[ED1,ED2,EI1,EI2]);

sinc(x):=sin(x)/x;
epsM:(k/q)**2*c**2+4*%pi*J0*d/(%i*q*(d1*E01+d2*E02+
    ED1*d1*sinc((k1-k)*d1/2)*exp(%i*(k1-k)*d1/2)+
    EI1*d1*sinc((k1+k)*d1/2)*exp(-%i*(k1+k)*d1/2)+
    ED2*d2*sinc((k2-k)*d2/2)*exp(%i*(k2-k)*d2/2)+
    EI2*d2*sinc((k2+k)*d2/2)*exp(-%i*(k2+k)*d2/2)));

epsM(z):=subst(k=z,epsM);
epsMD2:at(diff(epsM(z),z,2),z=0.001);

mysqrt(x):= if imagpart(sqrt(x))<0 then -sqrt(x) else sqrt(x);
#mysqrt(x):= sqrt(x);
epsM:subst([k1=mysqrt(e1)*q/c,k2=mysqrt(e2)*q/c,d1=80.0,d2=20.0,d=100.0,J0=1.0],epsM(0.001));

eps2:1.0-1.0/(q**2+%i*q*0.01);
eps1:(1.645+5.0732*10**3*q**2/(2*%pi*c)**2 + 2.31*10**7*q**4/(2*%pi*c)**4)**2;

epsM:subst([e1=eps1,e2=eps2,c=197.33],epsM);

epsM1(x):=subst(q=x,epsM);
epsM2re(x):=float(realpart(epsM1(x)));
epsM2im(x):=float(imagpart(epsM1(x)));

epsMlist:makelist([i,epsM2re(i),epsM2im(i)],i,1.5,3.5,0.04);
epsOut:openw("epsMwithoutInt.dat");
write_data(epsMlist,epsOut);
close(epsOut);

[Gunter Königsmann]

I once had a calculation that if done with floats suffered from the problem that adding a small float to a big one doesn't change the big float. Redoing the calculation with exact numbers exceeded the amount of memory gcl can allocate. But a maxima that was compiled with clips did the trick - even if it took about 8 hours.

[Wolfgang Dautermann]

If I remember correctly, the Maxima startscript ($INSTALLPREFIX/bin/maxima or src/maxima.in in the Maxima sourcecode) limits GCL memory consumption; see commit https://sourceforge.net/p/maxima/code/ci/9e712a93

You can try to increase/remove this limit.

[Guillermo P. Ortiz]

Thanks to @dauti, I am noticing that environment variable can manage the memory used. By default the maxima setup this variable more aggressive in order to allows multiple task. May be I can probe to increase this values, but not sure.
I have by default
export GCL_GC_PAGE_THRESH=0.2
export GCL_GC_ALLOC_MIN=0.01
export GCL_GC_PAGE_MAX=0.05
and
export GCL_MEM_MULTIPLE=0.2

what values for each would you suggest me?
Regards

[Robert Dodier]

Hi Guillermo, it would extremely helpful to help you debug what's going on if you post the data file, or a small part of it (with an indication of how big the whole file is), or a made-up file which has the same format. I agree with the general observations about numeric vs symbolic, etc., but the details about the problem, and the path to getting to a solution, depend on the data.

[Guillermo P. Ortiz]

Hi Robert,
I understand you. Here, for simplicity, consider the error output

+++++++
Maxima encountered a Lisp error:

Condition in MACSYMA-TOP-LEVEL [or a callee]: INTERNAL-SIMPLE-ERROR: The storage for CONS is exhausted. 1614811 pages allocated. Use ALLOCATE to expand the space.

Automatically continuing.
++++++

obtained for the simpler case code that i post here

https://sourceforge.net/p/maxima/bugs/3933/#3743

without intepolated data file issue. As you can see, it is not a misleading on how to
manage the interpolation

Regards

[Guillermo P. Ortiz]

Sorry,
here the bug report related to

https://sourceforge.net/p/maxima/bugs/3933/#15dd

Please report bugs to:
https://sourceforge.net/p/maxima/bugs
To report a bug, you must have a Sourceforge account.
Please include the following information with your bug report:

---

Maxima version: "5.44.0"
Maxima build date: "2021-04-24 14:52:58"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.12"
User dir: "/home/gortiz/.maxima"
Temp dir: "/tmp"
Object dir: "/home/gortiz/.maxima/binary/5_44_0/gcl/GCL_2_6_12"
Frontend: false

---

The above information is also reported by the function 'build_info()'.

[Robert Dodier]

Guillermo, I've worked with the script which you posted above (https://sourceforge.net/p/maxima/bugs/3933/#3743). I made a few changes and it seems to be working as intended now. I'll attach the modified script, a patch to show just what changed, and the output I get.

I made two kinds of changes: (1) postpone calling mysqrt and sinc until they have numerical arguments. There are different ways to achieve that; I just quoted the function calls, i.e. writing 'mysqrt(...) instead of mysqrt(...).

(2) Apply rectform to reduce numerical expressions involving complex numbers to just nnn + %i*nnn. This is kind of tricky -- I had to try putting rectform in different places in order to get a simple result. For a few attempts, I forget exactly what I did, but anyway the results were large, messy expressions.

[Guillermo P. Ortiz]

Thank you very much for your patch, Robert.
Just I test your version for mysqrt() and it seems working fine.
But I am not sure i rectform are using the default sqrt instead of mysqrt, it could be?
The question arise from the fact that I am not specting a negative sign of the imaginary part of the output as one obtained here https://sourceforge.net/p/maxima/bugs/_discuss/thread/e6e3db30f2/a456/attachment/dodier-epsMwithoutInt.dat

Regards

[Robert Dodier]

Guillermo, I've worked on the script some more, and I am getting some different values for the output. Do you want to take a look at the output file dodier-epsMwithoutInt.dat-take-2 ?

I modified the script some more to try harder to ensure that the sqrt is the principal argument sqrt. The whole business is pretty complex, and pretty fragile, so I'm not convinced it is working right. It would be a lot easier if Maxima had an option to give the principal argument. Anyway I would be interested to hear if the results are getting any closer to what you were expecting.

[Guillermo P. Ortiz]

Hello Robert,
I appreciate your effort very much, but unfortunately this result has
unexpected results also. The imaginary part should be positive.
Let me ask for some guide to understand maxima syntax used in your work
1) What is the meaning of the first two lines? linel:100 and
display2d:false
2) what is %emode and why is located there first as false and after as true
3)mathdeclare I suspect that is for aa. but you are using a lot of other
function like lambda, constantp and apply1
I guess that you are trying to take a countercloclwise in defining
exp(i*phi) for all phi, also when phi<0, is not it?

Regards

[Stavros Macrakis]

Bug reports are not the best place to ask questions about Maxima in general.
For reference information about Maxima functions, please read the documentation. Within Maxima, you can enter, for example,

? linel
? matchdeclare

etc.

[Robert Dodier]

Guillermo, the %emode flag and stuff about applying the rule to enforce counterclockwise arguments was an attempt to subvert Maxima's tendency to introduce negative (clockwise) arguments; as you have observed, that was maybe not entirely successful.

I am still looking at this problem, since the problem of trying to work with only positive (counterclockwise) arguments seems more general. I'll keep you posted.