Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I tried to make a workable example here. The problem is that it is very very slow. I must be doing something very stupid out here. Can somebody give me some clue how to evaluate this multiple sum faster. Right now it takes around 2 mins. However, for real problem I need the parameter Npol to be of the order of ~10.

Clear["Global`*"];
SetSystemOptions["CompileOptions"->"CompileReportExternal"->False];
DRSH=N[(94688879/10000000)*10^(-7)];
Npol=2;
nk=40;
PN=Evaluate[Table[0.0025 n,{n,1,Npol}]];
OmN=Evaluate[Table[(1.0/DRSH)^((-n+2)/Npol),{n,1,Npol}]];
fpol[n_, On_]:=(-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n)/; On== OmN[[1]]
fpol[n_, On_]:=(-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n) (0.77 2.0^-n-On)/(0.77 2.0^-n+On)

OverLaps1=Compile[{{k, _Integer},{n,_Real}, {m,_Real}},Exp[-0.5 Sum[If[q!= k,(fpol[q, n]-fpol[q, m])^2,0],{q,0,nk}]],CompilationTarget->"C"];
OverLaps2=Compile[{{k, _Integer},{n,_Real}, {m,_Real}},Exp[-0.5 Sum[If[q!= k,(fpol[q, n]+fpol[q, m])^2,0],{q,0,nk}]],CompilationTarget->"C"];
NpolSum=Compile[{{i, _Integer},{j, _Integer}, {k, _Integer}}, (Sum[ Sum[PN[[n]]  PN[[m]] fpol[k, OmN[[n]]]^i fpol[k, OmN[[m]]]^j Exp[-0.5 (fpol[k, OmN[[n]]]^2+ fpol[k, OmN[[m]]]^2)] ( OverLaps1[k, OmN[[n]], OmN[[m]]]+ (-1.0)^j   OverLaps2[k, OmN[[n]], OmN[[m]]]),{n,1,Npol}],{m,1,Npol}])^2,CompilationTarget->"C"];

S2Entrp=Compile[{{k, _Integer}},Sum[Sum[ ( 1/(i! j!)  )  NpolSum[i, j, k],{i,0,10}],{j,0,10}]];
ListLogLinearPlot[Table[{k,S2Entrp[ k]}, {k,0, nk}],PlotStyle->Directive[PointSize[Medium],Red], PlotRange-> All]//AbsoluteTiming

I thank you for your time.

share|improve this question
    
Without looking at the code in detail, did you try replacing Sum[...,{}] by Total@Table[...,{}]? –  acl Mar 29 '13 at 18:59

1 Answer 1

up vote 1 down vote accepted

One significant thing you could do is memoize the definitions for fpol and both OverLaps functions, as these are called multiple times with the same arguments.

I would also suggest creating a definition directly for fpol[n_, OmN[[1]]] rather than using the condition on the RHS.

I have ditched the Compile here, as I don't think it's giving any benefit due to all the calls to MainEvaluate (BTW suppressing the warnings doesn't make the problem go away...)

fpol[n_, OmN[[1]]] := fpol[n, OmN[[1]]] = (-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n)
fpol[n_, On_] := fpol[n, On] = (-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n) (0.77 2.0^-n - On)/(0.77 2.0^-n + On)

OverLaps1[k_, n_, m_] := OverLaps1[k, n, m] = Exp[-0.5 Sum[If[q != k, (fpol[q, n] - fpol[q, m])^2, 0], {q, 0, nk}]]    
OverLaps2[k_, n_, m_] := OverLaps2[k, n, m] = Exp[-0.5 Sum[If[q != k, (fpol[q, n] + fpol[q, m])^2, 0], {q, 0, nk}]]

NpolSum[i_, j_, k_] := (Sum[Sum[PN[[n]] PN[[m]] fpol[k, OmN[[n]]]^i fpol[k, OmN[[m]]]^j 
Exp[-0.5 (fpol[k, OmN[[n]]]^2 + fpol[k, OmN[[m]]]^2)] (OverLaps1[k, OmN[[n]], 
         OmN[[m]]] + (-1.0)^j OverLaps2[k, OmN[[n]], OmN[[m]]]), {n, 1, Npol}], {m, 1, Npol}])^2

S2Entrp[k_] := Sum[Sum[(1/(i! j!)) NpolSum[i, j, k], {i, 0, 10}], {j, 0, 10}]

I think there are further improvements that could be made, for example making fpol and OverLaps Listable and replacing NpolSum with something nifty using Outer, but the memoization appears to be the key to increasing the performance here.

With these changes the plot takes under half a second to create:

ListLogLinearPlot[Table[{k, S2Entrp[k]}, {k, 0, nk}], 
 PlotStyle -> Directive[PointSize[Medium], Red], PlotRange -> All] // AbsoluteTiming

enter image description here

share|improve this answer
    
Yes,thank you. I also found out that Compile doesn't help here. Again thank you for the quick tip. –  SSB Mar 29 '13 at 23:43

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.