Generating Gaussian random fields

This Mathematica notebook accompanies the paper “Local, algebraic simplifications of Gaussian random fields” by Theodor Bjorkmo and M. C. David Marsh. It is structured as follows:

1. Function definitions and graph settings
2. Setup
3. Generation of random function and evaluation of (Subscript[E, αβ]) matrices
4. Example: hyperparameter estimation
5. Example: plotting d=2 potentials

The code is straightforward to use. Section 1 need not be touched, and in section 2, the user only needs to specify the dimension of the space (d), order of expansion ((Subscript[n, max])), and the number of parallel kernels to be used. Section 3 generates the random function (f), with a seed that can be changed, and evaluates the (Subscript[E, αβ]) matrices for the specific d and (Subscript[n, max]).

Section 4 gives an example of hyperparameter estimation for GRFs. By maximising the log likelihood, the code tries to find the values of l and h set in the beginning of the section, and plots the 68% and 95% probability contour lines.

Section 5 is a short cell that plots the function in the case d=2.

TB is funded by an STFC studentship at DAMTP, University of Cambridge. DM is funded by Stephen Hawking Advanced Fellowship from the Centre for Theoretical Cosmology, DAMTP, University of Cambridge.

Function definitions and graph settings

(*Evaluates the matrix elements Dαβ*)
DαβFn[in_]:=Times@@(Factorial/@Tally[in][[All,2]]);

EαβPermListFn[baselist_,list_]:=DeleteDuplicates[Sort[#,Greater]&/@Flatten[Outer[Join,list,baselist,1,1],1]];

(*Evaluates the matrix elements Eαβ*)
EαβFn[in1_,in2_]:=Block[{int,count, mincount,overlap,c,sgn,remainder,tally1,t1,t2,tab1,tab2,tab},

t1=Tally[in1];
t2=Tally[in2];
tab1=Table[0,{i,d}];
tab2=Table[0,{i,d}];
Do[tab1[[t1[[i,1]]]]=t1[[i,2]],{i,Length[t1]}];
Do[tab2[[t2[[i,1]]]]=t2[[i,2]],{i,Length[t2]}];
tab=tab1-tab2;

If[AllTrue[tab,EvenQ[#]&&#≥0&],
c=Fold[#1Binomial[#2[[1]],#2[[2]]](#2[[1]]-#2[[2]]-1)!!&,1,Thread[{tab1,tab2}]];
If[OddQ[(Length[in1]-Length[in2])/2],
sgn=-1;,sgn=1;];
Return[sgn c]
,Return[0]
];
];


(*The two below functions evaluate the value of the random function f without expanding the Taylor coefficients as tensors*)
contrFn[ct_,k_,Φ_]:=Block[{t1=ct,vec=Φ,coeffOrder=k,out,iListInt,tally,ptable},
out=0.0;
iListInt=iList[[coeffOrder+1]];
Do[
tally=Tally[iListInt[[i]]];
ptable=Table[((vec[[tally[[j,1]]]])^(tally[[j,2]])),{j,1,Length[tally]}];

out=out+t1[[i]]*Multinomial[Sequence@@(tally[[1;;-1,2]])]*Product[ptable[[j]],{j,Length[ptable]}];
,{i,1,Length[iListInt]}];
out
];

f[Φ_]:=Block[{vec=Φ,out},
out=0.0;
Do[
out=out+contrFn[fα[[i+1]],i,vec]/Factorial[i]//N
,{i,0,order}];
out];

Needs["MaTeX`"]
Labels2=Directive[Black,FontFamily->"Latin Modern Roman",FontSize->12];
Colors3=Thread[List[{,,,,,,,,,,,,,,,,,,,},Thick]];
SetOptions[MaTeX,(*"Preamble"->{"\\usepackage{color,txfonts}"},*)FontSize->16]
hexList={"#A15562","#A45C73","#A36685","#9E7096","#957CA6","#8888B2"}; (*Purple/blue*)
hexToRGB=RGBColor@@(IntegerDigits[#~StringDrop~1~FromDigits~16,256,3]/255.)&;

Setup (choose dimension, expansion order, and number of parallel kernels)

(*Dimension of the GRF*)
d=4;
(*Highest order in the Taylor expansion*)
order=5;

kernels=32;






(*Number of coefficients at each order*)
di[x_]:=(d+x-1)!/((d-1)!x!);

dList=Map[di,Range[0,order]];

(*dList[[6]]Total[dList[[7;;]]]*)

(*Generally, e = even and o = odd*)
edList=dList[[1;;;;2]];
odList=dList[[2;;;;2]];
edListAcc=Accumulate[edList];
odListAcc=Accumulate[odList];

enList=Range[0,order][[1;;;;2]];
onList=Range[0,order][[2;;;;2]];


(*Creating the list of unique coefficients*)
iListFn[baselist_,list_]:=DeleteDuplicates[Sort[#,Greater]&/@Flatten[Outer[Join,list,baselist,1,1],1]];
iList=FoldList[iListFn[#2,#1]&,Table[{i},{i,d}],Table[Table[{i},{i,d}],{i,order-1}]];

eList=iList[[2;;;;2]];
eList=Join[{{{}}},eList];
oList=iList[[1;;;;2]];

iList=Join[{{{}}},iList];

eIndexList=Flatten[eList,1];
eIndexList[[1]]={{}};
oIndexList=Flatten[oList,1];

eOrder=Floor[order,2]/2+1;
oOrder=order+1-eOrder;

Do[
Do[
iListPosFn[iList[[i,j]]]=j;
,{j,1,Length[iList[[i]]]}];
,{i,1,order+1}];

Generate random function f and compute the (Subscript[E, αβ]) matrices

(*Random seed that can be changed*)
seed=3;


CloseKernels[];
LaunchKernels[kernels];

SeedRandom[seed];

Print["Even derivatives start at "<>DateString[]];

AbsoluteTiming[
eEαβList={};
eμαList={};
eCoeffs=Length[eIndexList];
eμα=ConstantArray[0,eCoeffs];
eMonitor=0;
Monitor[
efα=Reap[Do[
(*Print@i;*)
Block[{i1,i2,μαint,σList,fαint,Eαβint,EαβVecint,l1,l2,perms,vals,positions,posList,valList,ruleList},
eMonitor=i;
μαint=eμα[[1;;edList[[i]]]];

eμα=Drop[eμα,edList[[i]]];
σList=Sqrt/@DαβFn/@eList[[i]];

Sow[
fαint=Table[RandomVariate[NormalDistribution[μαint[[j]],σList[[j]]]],{j,1,Length[σList]}]];

If[i<eOrder,
l1=eList[[i]];
ruleList=Flatten[ParallelTable[Block[{steps,intList,iiTable,pos1},
steps=(2(eOrder-1)-Length[j])/2;
iiTable=Table[{i,i},{i,d}];
pos1=iListPosFn[j];
If[steps>0&&IntegerQ[steps],
intList={j};
Flatten[Reap[Do[
intList=EαβPermListFn[iiTable,intList];
vals=Map[EαβFn[#,j]&,intList];
positions=Map[{iListPosFn[#]+edListAcc[[i+k-1]]-edListAcc[[i]],pos1}&,intList];
Sow[Thread[positions->vals]],{k,steps}]][[2,1]],1]
]],{j,l1},Method->"CoarsestGrained"],1];

(*If[i==2,Print[ruleList]];*)
Eαβint=SparseArray[ruleList,{edListAcc[[eOrder]]-edListAcc[[i]],edList[[i]]}];


(*If[i==2,Print[Normal[Eαβint]]];*)
eμα=eμα+Eαβint.(fαint-μαint);
];





AppendTo[eμαList,μαint];
If[i<eOrder,
AppendTo[eEαβList,Eαβint];
];
]
,{i,1,eOrder}]][[2,1]];

,Row[{"eOrder: "<>ToString[eMonitor]<>" of "<>ToString[eOrder]}]];


Print["Odd derivatives start at "<>DateString[]];

(*SeedRandom[seed];*)

oMonitor=0;
oEαβList={};
oμαList={};
oCoeffs=Length[oIndexList];
oμα=ConstantArray[0,oCoeffs];
Monitor[
ofα=Reap[Do[
Block[{i1,i2,μαint,σList,fαint,Eαβint,l1,l2,perms,vals,positions,posList,valList,ruleList},
(*Print@i;*)
oMonitor=i;

μαint=oμα[[1;;odList[[i]]]];

oμα=Drop[oμα,odList[[i]]];
σList=Sqrt/@DαβFn/@oList[[i]];
Sow[
fαint=Table[RandomVariate[NormalDistribution[μαint[[j]],σList[[j]]]],{j,1,Length[σList]}]
];


If[i<oOrder,
l1=oList[[i]];

ruleList=Flatten[ParallelTable[Block[{steps,intList,iiTable,pos1},
steps=(2oOrder-1-Length[j])/2;
iiTable=Table[{i,i},{i,d}];
pos1=iListPosFn[j];
If[steps>0&&IntegerQ[steps],
intList={j};
Flatten[Reap[Do[
intList=EαβPermListFn[iiTable,intList];
vals=Map[EαβFn[#,j]&,intList];
positions=Map[{iListPosFn[#]+odListAcc[[i+k-1]]-odListAcc[[i]],pos1}&,intList];
Sow[Thread[positions->vals]],{k,steps}]][[2,1]],1]
]],{j,l1},Method->"CoarsestGrained"],1];

Eαβint=SparseArray[ruleList,{odListAcc[[oOrder]]-odListAcc[[i]],odList[[i]]}];


oμα=oμα+Eαβint.(fαint-μαint);
];




AppendTo[oμαList,μαint];

If[i<oOrder,
AppendTo[oEαβList,Eαβint];
];
]
,{i,1,oOrder}]][[2,1]];
,Row[{"oOrder: "<>ToString[oMonitor]<>" of "<>ToString[oOrder]}]];

][[1]]



fα=Riffle[efα,ofα];
μαList=Riffle[eμαList,oμαList];



fαProj=fα[[All,1]];

fProj[x_]:=Sum[(1/(i!))fαProj[[i+1]](x^i),{i,0,order}];

limit=5;

Plot[fProj[x],{x,-limit,limit},PlotRange->{{-limit,limit},{-5,5}}]

CloseKernels[];


eEαβListFn[i_,j_]:=eEαβList[[j,edListAcc[[i-1]]+1-edListAcc[[j]];;edListAcc[[i]]-edListAcc[[j]]]];
oEαβListFn[i_,j_]:=oEαβList[[j,odListAcc[[i-1]]+1-odListAcc[[j]];;odListAcc[[i]]-odListAcc[[j]]]];
eDαβList=Table[DiagonalMatrix@(DαβFn/@eList[[i]]),{i,eOrder}];
oDαβList=Table[DiagonalMatrix@(DαβFn/@oList[[i]]),{i,oOrder}];
eDαβInvList=Table[DiagonalMatrix@(((DαβFn[#])^(-1))&/@eList[[i]]),{i,eOrder}];
oDαβInvList=Table[DiagonalMatrix@(((DαβFn[#])^(-1))&/@oList[[i]]),{i,oOrder}];

Example: hyperparameter estimation

(*Setting the values for l and h. g is the rescaled version of f*)
lVal=0.33;
hVal=1.22;
gα=fα hVal Table[(lVal^(-i)),{i,0,order}];
egα=gα[[1;;;;2]];
ogα=gα[[2;;;;2]];

(*This code is written with fixed m=0, and computes the log likelihood in a region of hyperparater space set below. The isoprobability lines at 68% and 95% are plotted at the end of this cell. *)

(*Setting the grid over which the likelihood will be evaluated*)
lEstMin=0.25; 
lEstMax=0.4; 
hEstMin=0.05; 
hEstMax=3; 
gridSteps=100;

mEst=0;


CloseKernels[];
LaunchKernels[kernels];

grid=ParallelTable[Block[{lInt,hInt,eμαIntList,oμαIntList,eΔαIntList,oΔαIntList,lnP},
lInt=lEstMin+((i-1)/(gridSteps-1))(lEstMax-lEstMin);
hInt=hEstMin+((j-1)/(gridSteps-1))(hEstMax-hEstMin);
eμαIntList={{mEst}};
oμαIntList={Table[0,d]};
Do[Block[{int},
int=Total@Table[(lInt^(2(k-l)))eEαβListFn[l,k].(egα[[k]]-eμαIntList[[k]]),{k,1,l-1}];
AppendTo[eμαIntList,int];];
,{l,2,eOrder}];
Do[Block[{int},
int=Total@Table[(lInt^(2(k-l)))oEαβListFn[l,k].(ogα[[k]]-oμαIntList[[k]]),{k,1,l-1}];
AppendTo[oμαIntList,int];];
,{l,2,oOrder}];
eΔαIntList=egα-eμαIntList;
oΔαIntList=ogα-oμαIntList;

lnP=Block[{even,odd},
even=Total@Table[-(1/2)(hInt^(-2))(lInt^(2enList[[i]]))eΔαIntList[[i]].eDαβInvList[[i]].eΔαIntList[[i]]-edList[[i]](Log[hInt]-enList[[i]]Log[lInt])-((edList[[i]])/2)Log[2π]-(1/2)Tr[Log[eDαβList[[i]]]],{i,1,eOrder}];
odd=Total@Table[-(1/2)(hInt^(-2))(lInt^(2onList[[i]]))oΔαIntList[[i]].oDαβInvList[[i]].oΔαIntList[[i]]-odList[[i]](Log[hInt]-onList[[i]]Log[lInt])-((odList[[i]])/2)Log[2π]-(1/2)Tr[Log[oDαβList[[i]]]],{i,1,oOrder}];
even+odd];

{lInt,hInt,lnP}

],{i,gridSteps},{j,gridSteps},Method->"CoarsestGrained"];
CloseKernels[];
gridMax=MaximalBy[Flatten[grid,1],#[[3]]&][[1]];
Print["Correct values: ", {lVal,hVal}];
Print["Best guess: ",gridMax[[1;;2]]];

grid2=Flatten[grid,1];
latticepoints=Length[grid2];
p0grid={#[[1]],#[[2]],Exp[#[[3]]]}&/@grid2;
ptot= Total[p0grid[[;;,3]]];
pgrid={#[[1]],#[[2]],#[[3]]/ptot}&/@p0grid;


pSorted =Sort[pgrid,#1[[3]]>#2[[3]]&];
pAccumulated=Accumulate[pSorted[[;;,3]]];
P68pos=FirstPosition[pAccumulated,x_/;x≥ 0.68][[1]]-1;
P95pos=FirstPosition[pAccumulated,x_/;x≥ 0.95][[1]]-1;
P68iso=pSorted[[P68pos,3]];
P95iso=pSorted[[P95pos,3]];

n1=3;
plot1=ListContourPlot[pSorted, PlotRange->Full(*{{lEstMin,lEstMax},{hEstMin,hEstMax}}*),(*Contours->{0.0013},*)FrameLabel->MaTeX/@{"l","h"},Contours->{{P95iso,{Colors3[[n1]],Dashed}},{P68iso,Colors3[[n1]]}},ContourShading->None,LabelStyle->Labels2]

Example: plotting 2d potentials

CloseKernels[];
LaunchKernels[kernels];

limit=4.0;

points=50;

If[d==2,
coordList=Table[(((2i)/points) -1)limit,{i,0,points}];

AbsoluteTiming[
plot3dList=Flatten[Parallelize[Outer[{#1,#2,f[{#1,#2}]}&,coordList,coordList],Method->"CoarsestGrained"],1];]//Print;



landscapePlot=ListPlot3D[plot3dList,ClippingStyle->None,LabelStyle->None,TicksStyle->None,ColorFunction->Function[{x,y,z},Blend[hexToRGB/@hexList,z]],Axes->True,Boxed->False,AxesOrigin->{0,0,0},AxesStyle->Opacity[0],TicksStyle->None];
landscapePlot//Print;

];
CloseKernels[];