In [8]:
import cft
from numpy import random, fft, linalg, pi
import numpy as np
from IPython.display import display, Image, HTML
from gnuplot import plot, plot_three, plot_six, plot_pdf, Multiplot
N = 256
L = 22
B = cft.Box(2, N, L)
Constrained Gaussian Random fields can be as easy as pie. What do we know about GRFs? We have a one-point function of
$$\mathcal{P}(f_1) = \frac{1}{\sqrt{\tau v}} \exp \left[-\frac{f_1}{2v}\right]$$and an $N$-point function of
$$\mathcal{P}(f) = \frac{1}{\sqrt{\tau^N \det M}} \exp\left[-\frac{1}{2}f^{\dagger}M^{-1}f\right],$$where $f$ is now an $N$-vector of values, and $M$ the $(N \times N)$-covariance-matrix,
$$M_{ij} = \langle f_i^*f_j\rangle.$$In the case of Gaussian random fields, all the $N$-point functions can be described in terms of the two-point function, and we can write the covariance matrix as
$$M_{(2)} = \begin{pmatrix} v & \xi \\ \xi & v \end{pmatrix},$$where $\xi$ is the two-point correlation function. In the paper by van de Weygaert & Bertschinger (1996) (WB96), the $N$-point function is written in a functional form. We will be focussing on the computation of CGRFs, we will keep the matrix the notation for finite $N$. Computing the expectation value of a quantity $A(f)$,
$$\langle A\rangle = \int A(f) \mathcal{P}(f) {\rm d}^N f,$$assuming $\mathcal{P}(f)$ is properly normalised.
To show the difference between an uncorrelated and correlated random field, I have created three instances of a GRF below with power spectrum of $P(k) = k^{n}$, where $n = 0, -1, -2$. The value in each cell is an element in the vector $f$.
In [2]:
def garfield(B, wn, P):
f = fft.ifftn(fft.fftn(wn) * np.sqrt(P(B.K))).real
f /= f.std()
return f
P0 = cft.Power_law(2)
P1 = cft.Power_law(-1) * cft.Scale(B, 0.1)
P2 = cft.Power_law(-2)
random.seed(17)
wn = random.normal(0, 1, [N, N])
plot_three(B, wn, garfield(B, wn, P1), garfield(B, wn, P2))
We will influence the generating of a random field by imposing constraints. For example, we might want a density peak in the center of our box, then we can encode this wish in a series of constraints: scaled with a certain gaussian filter the gradient is zero, and the density has some positive value, while the second derivatives are all negative. The field is then subject to a set of $M$ constraints such that,
$$\Gamma = \left\{ C_i(f) = g_i ; i = 1, \dots, M \right\}.$$For practical purposes we have to make sure that each costraint function $C_i$ is linear in the sense that it can be expressed in terms of a linear combination of all elements in $f,$
$$C_i(f) = \langle\mathcal{C}_i, f\rangle,$$or if you will, there is a matrix $N \times M$ matrix $C$ that transforms $f$, to the set of constraint values $g$.
In particular the case where the constraint can be expressed as a convolution is common
$$C_i(f, x) = \frac{1}{N}\sum_i g(x - y_i) f(y_i).$$The problem is how to sample the possible constraint realisations properly,
$$\mathcal{P}\big(f|\Gamma\big) = \frac{\mathcal{P}\big(f \cap \Gamma\big)}{\mathcal{P} \big(\Gamma\big)} = \frac{\mathcal{P}\big(f\big)}{\mathcal{P}\big(\Gamma\big)}.$$Since the coefficients $c_i$ are linear combinations of Gaussian variables, they are themselves distributed as a multivariate Gaussian with the covariance matrix $Q_{ij} = \langle g_i^* g_j\rangle = CMC^{\dagger},$ (how do we show this last equality? $(AB)^{-1} = B^{-1}A^{-1}$ then substitute in $N$-point distribution)
$$\mathcal{P}\big(\Gamma\big) = \frac{1}{\sqrt{\tau^M \det Q}} \exp\left[-\frac{1}{2}g^{\dagger}Q^{-1}g\right].$$$$\mathcal{P}\big(f|\Gamma\big) = \sqrt{\frac{\tau^M \det Q}{\tau^N \det M}} \exp\left[-\frac{1}{2}\left(f^{\dagger}M^{-1}f - g^{\dagger}Q^{-1}g\right)\right]$$The term in the exponential can be written in the form $1/2F^{\dagger} M^{-1} F$, where $F = f - \bar{f}$. Defining $\bar{f}$, the mean field under the constraints as
$$\bar{f} := MC^{\dagger}Q^{-1}g = \big\langle f | \Gamma \big\rangle.$$The combination $MC^{\dagger}$ equals the cross-correlation between the random field and the set of constraints $\langle f C_i(f)\rangle$. If we were to substitute $C^{-\dagger}M^{-1}C^{-1}$ into this definition, we would get the expression $C^{-1} g$, but $C$ is in general not invertible. According to Bertschinger (1987) (B87), $f = \bar{f}$ is a stationary point of the action: $\delta S/\delta f = 0$ for $f = \bar{f}$. Where the action $S$ is the term in the exponential of the distribution function. This computation is apparently explained in Bardeen et al. (1986) (BBKS).
Moving back to WB96 and its appendix C, we now see that the constrained field is the sum of an average part $\bar{f}$ and a residual part $F$. The insight needed is that an unconstrained field $\tilde{f}$ will have non-zero constraint coefficients $\tilde{g} = C\tilde{f}$. Since everything is linear, all we need to do is add to the unconstrained field an average field with coefficients $g - \tilde{g}$. Suppose there are two sets of constraints $\Gamma_1$ and $\Gamma_2$, and we have a realisation $f_1$ with constraint coefficients $g_1 = C f_1$, then we can transform this realisation into $f_2$ with $g_2 = C f_2$ as follows
$$f_2 = f_1 - \bar{f}_1 + \bar{f}_2.$$This implies that there is a bijection between the realisations of two sets of constraints. In other words, for each realisation $f$, drawn from the distribution with constraints $\Gamma_1$, there is one with equal probability from $\Gamma_2$, or
$$\mathcal{P}\big(f_1 | \Gamma_1\big) = \mathcal{P}\big(f_2 | \Gamma_2\big).$$This also means that the residual field $F$ is independent of the chosen constraint coefficients.
Computing the average field seems like a lot of work, since there is the matrix $M$ in the expression. However, if we do all our work in Fourier space, the different modes $\hat{f}(k)$ are independent, and $M$ is diagonal. Again, because the Fourier transformation is linear, we retain all nice properties that we just derived.
In [3]:
import gnuplot
import imp ; imp.reload(gnuplot)
sigma = 1.0
S = cft.Cutoff(B) * cft.Scale(B, sigma) * cft.Pos([L/2, L/2])
extra = [(a * S) / (a * S).abs(B, P1) for a in
[cft.Potential() * cft.D(1), cft.Potential() * cft.D(2)]]
H = cft.make_peak(B, P1, S) + extra
gnuplot.plot_six(B, *[fft.ifftn(h(B.K)).real for h in H][:6])
Now we take a normal unconstrained field and compute the resulting coefficients.
In [4]:
cgrf = cft.CGRF(B, P1, H)
cgrf.generate_noise()
cgrf.set_coefficients(np.r_[1.0,0,0, -0.53, -0.53,0, 0, 0]*-4)
A1 = cgrf.triplet()
#cgrf.generate_noise(seed=53)
A2 = cgrf.triplet(cft.Potential())
plot_three(B, *A1)
print(cgrf.p_value(cgrf.g))
In [5]:
from scipy.interpolate import RectBivariateSpline
from functools import reduce
x = np.mgrid[-B.L/2:B.L/2:B.N*1j]
Q = np.mgrid[-B.L/2:B.L/2:128j,-B.L/2:B.L/2:128j]
if False:
for i in range(3):
iPhi = RectBivariateSpline(x, x, A2[i])
sphi = iPhi.ev(Q[0].flat,Q[1].flat).reshape([128,128])
X = [(Q - np.gradient(D * sphi, B.L/128.)).transpose([1,2,0]).reshape([128**2, 2]) \
for D in [0.2, 0.5, 1.5]]
cmd = ("set term pngcairo size 720,240 font 'Bitstream Charter, 9'", "set multiplot layout 1,3",
"unset key", "set xrange [-10:10] ; set yrange [-10:10]") + \
tuple(reduce(lambda x,y: x+y, (["plot '-' w dots", x] for x in X))) + \
("unset multiplot",)
plot(*cmd)
Q = np.indices(B.shape) * B.res
for i in range(3):
vx = fft.ifftn(fft.fftn(A2[i]) * -1j * B.K[0]).real
vy = fft.ifftn(fft.fftn(A2[i]) * -1j * B.K[1]).real
X = (Q + np.array([vx, vy])).transpose([1,2,0]) + L
if (i == 1):
X += random.normal(0, 1e-5, X.shape)
T0 = X.reshape([B.N**2,2])
T1 = np.roll(X, -1, axis=0); T1[-1,:,0] += L ; T1 = T1.reshape([B.N**2,2])
T2 = np.roll(X, -1, axis=1); T2[:,-1,1] += L ; T2 = T2.reshape([B.N**2,2])
T3 = np.roll(np.roll(X, -1, axis=0), -1, axis=1)
T3[-1,:,0] += L ; T3[:,-1,1] += L ; T3 = T3.reshape([B.N**2, 2])
tt1 = np.c_[T0, T1, T2]
tt2 = np.c_[T3, T2, T1]
f = open('tmp{0}.dat'.format(i), 'wb')
f.write("{0}\n".format(2*len(tt1)).encode())
np.savetxt(f, tt1)
np.savetxt(f, tt2)
f.close()
In [6]:
import os
f = open('tmp.dat', 'wb')
for i in range(3):
os.system("../nerve2/nerve nerve -i tmp{0} -N 512 -L 22 --txt tmp{0}.dat".format(i))
np.savetxt(f, A1[i].T)
f.write("\n\n".encode())
f.close()
In [10]:
M = Multiplot(3, 3, 2, 3)
plot_pdf("void-zeld",
M.set_term(),
"set xrange [-{0}:{0}] ; set yrange [-{0}:{0}]".format(L/2),
"set tics front",
"unset key", "unset colorbox", #"set size square",
"set multiplot",
"set cbrange [-2.5:2.5]",
"set xtics in format ''",
M.subplot(0, 0),
"plot 'tmp.dat' i {2} matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(N/L, L/2, 0),
M.subplot(0, 1),
"plot 'tmp.dat' i {2} matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(N/L, L/2, 1),
"set xtics format '%g'",
M.subplot(0, 2),
"set colorbox horizontal user origin {0},{1}-{3} size {2}, {3}".format(M.left(0) + M.ew*0.05, M.bottom(2) - M.eh/10, M.ew*0.9, M.eh/20),
"plot 'tmp.dat' i {2} matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(N/L, L/2, 2),
"unset colorbox",
"set cbrange [0.065:6]", "set log cb",
"set ytics in format ''", "set xtics in format ''",
M.subplot(1, 0),
"plot 'tmp{2}.nerve.conan' i 0 matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(512/L, L/2, 0),
M.subplot(1, 1),
"plot 'tmp{2}.nerve.conan' i 0 matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(512/L, L/2, 1),
M.subplot(1, 2),
"set xtics format '%g",
"set colorbox horizontal user origin {0},{1}-{3} size {2}, {3}".format(M.left(1) + M.ew*0.05, M.bottom(2) - M.eh/10, M.ew*0.9, M.eh/20),
"plot 'tmp{2}.nerve.conan' i 0 matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(512/L, L/2, 2),
"unset multiplot")
In [11]:
from IPython.display import Latex
def format_matrix(M):
s = "\\[\\left(\\begin{array}{" + "r"*M.shape[0] + "}\n "
s += " \\\\\n ".join([" & ".join(["{0: >6,.2f}".format(q) for q in row]) for row in np.array(M)])
s += "\n\\end{array}\\right)\\]\n"
return Latex(s)
format_matrix(cgrf.Q)
Out[11]:
In [12]:
np.sqrt(1/5)
Out[12]:
In [13]:
plot_three(B, A1[0], A1[2], A1[2] - A1[0])
In [14]:
cgrf.action(cgrf.g)
Out[14]:
In [ ]: