I am working with noncommuting objects and basically I am using Mathematica to sort large expressions by normal ordering and simplifying them. To this end, I first create an object called Boson
(These satisfy $[a,a^\dagger]=1$):
Clear[Boson, BosonC, BosonA]
Boson /: MakeBoxes[Boson[cr : (True | False), sym_], fmt_] :=
With[{sbox = If[StringQ[sym], sym, ToBoxes[sym]]},
With[{abox =
If[cr, SuperscriptBox[#, FormBox["\[Dagger]", Bold]], #] &@sbox},
InterpretationBox[abox, Boson[cr, sym]]]]
BosonA[sym_: String "a"] := Boson[False, sym]
BosonC[sym_: String "a"] := Boson[True, sym]
Next I alias the noncommutative product with CenterDot
as follows:
Unprotect[NonCommutativeMultiply];
Clear[NonCommutativeMultiply, CenterDot];
CenterDot[a__] := NonCommutativeMultiply[a];
NonCommutativeMultiply /:
MakeBoxes[NonCommutativeMultiply[a__], fmt_] :=
With[{cbox = ToBoxes[HoldForm[CenterDot[a]]]},
InterpretationBox[cbox, NonCommutativeMultiply[a]]]
Protect[NonCommutativeMultiply];
Clear[CRule]
CRule = {NonCommutativeMultiply[a_] :> a};
Then I define the Function
clean
which does the normal ordering. It is defined in terms of expand
. The main definition for expand
is
ClearAll@expand
SetAttributes[expand, HoldAll]
Unevaluated[
expand[expr_] :=
Block[
{NonCommutativeMultiply (*or times*)},
expr //. (*ReplaceRepeated instead of ReplaceAll*)
{times[left___, cnum_ /; FreeQ[cnum, (_Boson)],
right___] :> cnum*times[left, right],
times[left___,
cnum_ /; (! FreeQ[cnum, Times[n___?NumericQ, ___Boson]]),
right___] :>
Times @@ Apply[Power, Drop[FactorList[cnum], -1], 2]*
times[left, First[Last[FactorList[cnum]]], right],
times[left___, Boson[False, s_], Boson[True, s_], right___] :>
times[left, right] +
times[left, Boson[True, s], Boson[False, s], right],
times[left___, fst : Boson[_, s_], sec : Boson[_, t_],
right___] :>
times[left, sec, fst, right] /;
FreeQ[Ordering[{s, t}], {1, 2}], times[b_Boson] :> b,
times[] -> 1
}
]
] /. {HoldPattern[times] -> NonCommutativeMultiply }
More definitions for expand are
expand[Alternatives[NonCommutativeMultiply, CenterDot][
a1_, (a2_ + a3_)]] := expand[a1 ** a2] + expand[a1 ** a3]
expand[Alternatives[NonCommutativeMultiply, CenterDot][(a1_ + a2_),
a3_]] := expand[a1 ** a3] + expand[a2 ** a3]
Finally, the definition of clean
is
Clear[clean]
clean = Simplify[FixedPoint[expand, Distribute //@ #] //. CRule] &;
This has been working very well for computing quartic terms but now the problem is that I am computing higher order terms which can involve products of up to 16 oscillators. When I use the clean[]
function on those terms, the computation goes okay for a while but then the kernel starts using pretty much all the memory and it just stops doing anything else and I have quit the kernel. I am assuming that the Module[]
function is not deleting the temporary variables and is leaking memory. I tried setting $HistoryLength=0
but that didn't help either. What else can I do prevent this?
Update: I think the problem is in the implementation of FixedPoint
in clean
. If I manually apply expand
and Distribute
and then clear the cache and keep doing it till the computation doesn't take long (I take that as an indication that not much is changing), then I can finish the computation. This is a "dirty" fix but it would be nice to know if there is a way I can implement clearing the cache in clean
along with FixedPoint
.
Update: The data in one of the computations became so large that I kept running out of memory for everything and I had to keep dividing the expression into smaller and smaller parts to the point where it became hard to manage. Could be just that expressions are too big and I just don't have enough memory? Also edited the last block of code to remove a stray piece of code.
Update: Here is the full working notebook in pastebin.
Module
to be changed toBlock
? $\endgroup$Module
, the memory seems to be around 300-400 MB for a minute and then quickly shoots up to 2.5GB or so and stalls. WithBlock
it stays around 1.3GB and then shoots up and stalls. $\endgroup$NonCommutativeMultiply
bytimes
in order to avoid issues with the attributeFlat
? To avoid this, you can simply useexpand[expr_] := Block[{NonCommutativeMultiply}, replacementCode]
, which temporarily disregards any attributes (as any other definitions) ofNonCommutativeMultiply
. $\endgroup$