Initial commit

random_fields.py

+import os
+import itertools
+
+import numpy as np
+import scipy as sp
+
+from scipy.special import kv, gamma
+
+
+def eps():
+    return 1e-2
+
+
+def matern_cov(xs, ys, nu, sigma, l):
+    return matern_cov_dm(dist_matrix(xs, ys), nu, sigma, l)
+
+
+def matern_cov_dm(dist_mtx, nu, sigma, l):
+    s = (sigma**2) * np.power(2.0, 1 - nu) / gamma(nu)
+    ld = np.sqrt(2 * nu) * dist_mtx / l
+    kd = kv(nu, ld)
+
+    return s * np.power(ld, nu) * kd + eps() * np.eye(dist_mtx.shape[0])
+
+
+def dist_matrix(xs, ys):
+    n_pts = np.size(xs) * np.size(ys)
+    dist_mtx = eps() * np.ones((n_pts, n_pts), dtype=np.float32)
+
+    xys = list(itertools.product(xs, ys))
+    for row in range(len(xys)):
+        x, y = xys[row]
+        for col in range(row):
+            x1, y1 = xys[col]
+            dx = x - x1
+            dy = y - y1
+            dist = np.sqrt(dx * dx + dy * dy)
+            dist_mtx[row, col] = dist
+            dist_mtx[col, row] = dist_mtx[row, col]
+
+    return dist_mtx
+
+
+def make_matern_sampler(dist_mtx, n_sigma=eps()):
+    def sampler_impl(state, mean):
+        nu, sigma, l = state
+
+        cov_mtx = matern_cov_dm(dist_mtx, nu, sigma, l)
+
+        nc_inv = np.linalg.inv(cov_mtx + n_sigma * np.eye(cov_mtx.shape[0]))
+        cnc = np.matmul(cov_mtx, nc_inv)
+
+        mu = np.matmul(cnc, mean)
+        cov = cov_mtx - np.matmul(cnc, cov_mtx)
+        chol_cov = np.linalg.cholesky(cov)
+        n = np.random.normal(size=(cov_mtx.shape[0],))
+
+        return mu + np.matmul(chol_cov, n), cov
+
+    return sampler_impl
+
+
+def matern_parameter_posterior_nl(observed, field, field_cov, n_sigma=0):
+    field_diff = observed - field
+    nll = np.dot(field_diff, field_diff)
+
+    field_inv_cov = np.linalg.inv(field_cov)
+    field_nl_prior = np.matmul(np.matmul(field.T, field_inv_cov), field)
+
+    chol_cov = np.linalg.cholesky(field_cov)
+    tril = chol_cov[np.tril_indices(chol_cov.shape[0])]
+    log_det_prior = 2 * np.sum(tril)
+
+    return nll + field_nl_prior + 0.5 * (log_det_prior + chol_cov.shape[0] * np.log(2 * np.pi))
+
+
+def gibbs_metropolis(initial, proposal_sampler, f_nlp, field_sampler,
+                     n_samples=1000, burn_in=100, n_save=500):
+    chain = []
+    nlls = []
+    p_accs = []
+
+    state, observed_field = initial
+    field, field_cov = field_sampler(state, observed_field)
+    state_nll = f_nlp(state, field, field_cov)
+
+    for i in range(n_samples):
+        proposal_field, field_cov = field_sampler(state, observed_field)
+        proposal_state = proposal_sampler(state)
+        proposal_nll = f_nlp(proposal_state, proposal_field, field_cov)
+
+        acceptance = np.exp(state_nll - proposal_nll)
+
+        if np.random.random() < acceptance:
+            state = proposal_state
+            field = proposal_field
+            state_nll = proposal_nll
+            print("Gibbs-Metropolis [{}]: acceptance with p={}".format(i, acceptance))
+        else:
+            print("Gibbs-Metropolis [{}]: rejected with p={}".format(i, acceptance))
+
+        if i > (n_samples - n_save):
+            chain.append((state, field))
+
+        p_accs.append(acceptance)
+        nlls.append(state_nll)
+
+    return chain, p_accs, nlls
+
+
+
+
+
+
+
+
+

sample_matern.py

+import os
+import random
+import pickle
+
+import cv2
+import numpy as np
+
+import matplotlib
+from matplotlib import pyplot as plt
+
+import random_fields as rf
+
+
+def main():
+    w, h = 64, 64
+    xs = np.linspace(0, 5, w)
+    ys = np.linspace(0, 5, h)
+
+    n_pts = np.size(xs) * np.size(ys)
+
+    print("Pre-evaluating distance matrix")
+    dist_mtx = rf.dist_matrix(xs, ys)
+    print(dist_mtx.min())
+
+    dist_mtx_img = (dist_mtx - dist_mtx.min()) / (dist_mtx.max() - dist_mtx.min()) * 255
+    cv2.imwrite('./figs/dist_mtx.png', dist_mtx_img)
+
+    ns = [
+        np.random.normal(size=(n_pts,)),
+        np.random.normal(size=(n_pts,)),
+        np.random.normal(size=(n_pts,))
+    ]
+
+    print("Evaluating matern: nu=0.5, sigma=1, l=2")
+    k_smooth_05 = rf.matern_cov_dm(dist_mtx, nu=0.5, sigma=1, l=2)
+    print("Sampling...")
+    sample_and_save_fig(k_smooth_05, ns, w, h, 'k_smooth_05')
+
+    print("Evaluating matern: nu=1, sigma=1, l=2")
+    k_smooth_1 = rf.matern_cov_dm(dist_mtx, nu=1, sigma=1, l=2)
+    print("Sampling...")
+    field = sample_and_save_fig(k_smooth_1, ns, w, h, 'k_smooth_1')
+    with open('./data/field.pickle', 'wb') as f:
+        pickle.dump((xs, ys, field, 1, 1, 2), f)
+
+    print("Evaluating matern: nu=2, sigma=1, l=2")
+    k_smooth_2 = rf.matern_cov_dm(dist_mtx, nu=2, sigma=1, l=2)
+    print("Sampling...")
+    sample_and_save_fig(k_smooth_2, ns, w, h, 'k_smooth_2')
+
+
+def sample_and_save_fig(k_mtx, ns, w, h, prefix):
+    plt.figure()
+    plt.imshow(k_mtx)
+    plt.savefig('./figs/{}_cov_mtx'.format(prefix))
+    plt.close()
+
+    ck = np.linalg.cholesky(k_mtx)
+    field_sample = None
+    for i in range(len(ns)):
+        field_sample = np.matmul(ck, ns[i]).reshape((h, w))
+
+        plt.figure()
+        plt.imshow(field_sample)
+        plt.savefig('./figs/{}_{}'.format(prefix, i))
+        plt.close()
+
+    return field_sample
+
+
+if __name__ == '__main__':
+    random.seed(42)
+    np.random.seed(42)
+
+    main()
\ No newline at end of file