cloudy estimates added

assess_chain.py

+import os
+import argparse
+import pickle
+
+import numpy as np
+
+from scipy.special import kv, gamma
+from scipy.spatial import distance
+
+import matplotlib
+from matplotlib import pyplot as plt
+
+import random_fields as rf
+
+from external.multiESS import multiESS as m_ess
+
+
+def get_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--pickle_path", required=True,
+                        help="Path to pickled chain")
+    parser.add_argument("--outpath", required=True,
+                        help="Path where to save results")
+
+    return parser.parse_args()
+
+
+if __name__ == '__main__':
+    args = get_args()
+    w, h = 32, 32  # TODO load from pickles
+
+    with open(args.pickle_path, "rb") as f:
+        chain_res = pickle.load(f)
+
+        chain, p_accs, nlls = chain_res
+        theta_chain, field_chain = zip(*chain)
+
+        theta_chain = np.array(theta_chain)
+
+        mv_ess = m_ess.multiESS(theta_chain)
+        print("Multivariate ESS: {}".format(mv_ess))
+
+        with open(os.path.join(args.outpath, "mv_ess.txt"), "w") as mv_ess_f:
+            mv_ess_f.write("Multivariate ESS: {}".format(mv_ess))
+
+        figs_outpath = args.outpath
+        os.makedirs(figs_outpath, exist_ok=True)
+        rf.show_theta_chain(theta_chain, p_accs, nlls, figs_outpath, theta_content=["nu", "sigma", "l"])
+        field_mean, _ = rf.estimate_field_norm_posterior(field_chain, w, h, figs_outpath)
\ No newline at end of file

cloudy_estimate_matern.py

+import os
+import pickle
+import random
+import time
+
+import numpy as np
+
+import matplotlib
+from matplotlib import pyplot as plt
+
+import random_fields as rf
+
+
+def make_neg_log_posterior(noise_sigma, omxys, cloud_mask):
+    def impl(observed_field, field, state):
+        cloudy_field = field[cloud_mask.astype(np.bool)]
+        cloudy_observed = observed_field[cloud_mask.astype(np.bool)]
+
+        noise_cov = np.eye(np.size(cloudy_field)) * (noise_sigma ** 2)
+        nll = rf.neglog_mvnpdf(cloudy_field, cloudy_observed, noise_cov)
+
+        nu, sigma, l = np.exp(state)
+        cloudy_dist_mtx = rf.mdist(omxys, omxys)
+        field_cov = rf.matern_cov_dm(cloudy_dist_mtx, nu, sigma, l)
+        nl_prior = rf.neglog_mvnpdf(cloudy_field, 0, field_cov)
+
+        return nl_prior + nll
+    return impl
+
+
+def make_field_sampler(noise_sigma, omxys, cmxys, flat_cloud_mask):
+    def impl(observed_field, state):
+        theta = np.exp(state)
+
+        shrinked_observed_field = observed_field[flat_cloud_mask]
+        interpolated_field = rf.interpolate_field(shrinked_observed_field, theta, noise_sigma,
+                                                  omxys, cmxys)
+
+
+        return interpolated_field
+    return impl
+
+
+def make_cloud(w, h):
+    cloud_mask = np.ones((w, h), dtype=np.float32)
+
+    sx, ex = int(w / 2 - w / 8), int(h / 2 + h / 8)
+    sy, ey = int(w / 2 - w / 8), int(h / 2 + h / 8)
+    cloud_mask[sy:ey, sx:ex] = 0
+
+    return cloud_mask, sx, ex, sy, ey
+
+
+def main():
+    with open('./data/field.pickle', 'rb') as f:
+        xs, ys, field, nu, sigma, l = pickle.load(f)
+
+    w, h = np.size(xs), np.size(ys)
+
+    print("Pre-evaluating distance matrix")
+    dist_mtx = rf.dist_matrix(xs, ys)
+    flat_gt_field = np.ravel(field)
+
+    cloud_mask, sx, ex, sy, ey = make_cloud(w, h)
+    flat_cloud_mask = cloud_mask.ravel().astype(np.bool)
+
+    noise_sigma = 1e-2
+    noise = np.random.normal(scale=noise_sigma, size=flat_gt_field.shape)
+
+    observed_field = (flat_gt_field + noise) * flat_cloud_mask
+
+    omxs, omys = np.meshgrid(xs, ys)
+    # plt.figure()
+    # plt.imshow(omxs)
+    # plt.figure()
+    # plt.imshow(omys)
+    omxs = omxs.ravel()[flat_cloud_mask]
+    omys = omys.ravel()[flat_cloud_mask]
+    omxys = np.array([omxs, omys]).T
+
+    fmxs, fmys = np.meshgrid(xs, ys)
+    fmxs = fmxs.ravel()
+    fmys = fmys.ravel()
+    fmxys = np.array([fmxs, fmys]).T
+    field_sampler = make_field_sampler(noise_sigma, omxys, fmxys, flat_cloud_mask)
+
+    initial_theta = np.array([nu, sigma, l]) + 0.7
+    initial_state = np.log(initial_theta)
+
+    proposal_cov = 1e-4 * np.eye(np.size(initial_state))
+    lower = np.log(np.array([0.5, 0.5, 0.5]))
+    upper = np.log(np.array([3, 3, 3]))
+
+    neg_log_posterior = make_neg_log_posterior(noise_sigma, omxys, cloud_mask.ravel())
+
+    print("Starting MCMC")
+    chain_res = rf.gibbs_within_metropolis(
+        initial_state=initial_state,
+        initial_field=field_sampler(observed_field, initial_state),
+        observed_field=observed_field,
+        neg_log_posterior=neg_log_posterior,
+        field_sampler=field_sampler,
+        proposal_cov=proposal_cov,
+        lower=lower, upper=upper,
+        adaptation_start=1000
+    )
+
+    with open('./data/cloudy_matern_chain.pickle', 'wb') as f:
+        pickle.dump(chain_res, f)
+
+    chain, p_accs, nlls = chain_res
+    theta_chain, field_chain = zip(*chain)
+    figs_outpath = './figs/cloudy_matern/'
+    os.makedirs(figs_outpath, exist_ok=True)
+    rf.show_theta_chain(theta_chain, p_accs, nlls, figs_outpath, theta_content=["nu", "sigma", "l"])
+    field_mean, _ = rf.estimate_field_norm_posterior(field_chain, w, h, figs_outpath)
+
+    field_err = field - field_mean
+    plt.figure()
+    plt.imshow(field_err)
+    plt.colorbar()
+    plt.title("GT - mean estimate")
+    plt.savefig(os.path.join(figs_outpath, "field_err.png"))
+    plt.close()
+
+
+if __name__ == '__main__':
+    random.seed(42)
+    np.random.seed(42)
+
+    main()

docs/make_study_diary.sh

+#!/bin/bash
+
+./make_doc.sh "laine_study"
\ No newline at end of file

docs/texs/laine_study.tex

+\documentclass[a4paper,runningheads]{article}
+
+\usepackage{lineno,hyperref}
+
+\usepackage{graphicx}
+
+\usepackage{xcolor}
+\usepackage{amsmath, amssymb}
+\usepackage{graphicx}
+\usepackage{subfig}
+\usepackage[paperwidth=210mm,paperheight=297mm]{geometry}
+\usepackage{enumitem}
+\usepackage{mathtools}
+
+\renewcommand\UrlFont{\color{blue}\rmfamily}
+
+\DeclareMathOperator{\E}{\mathbb{E}}
+\DeclareMathOperator{\diag}{diag}
+\DeclareMathOperator{\sign}{sign}
+\DeclareMathOperator*{\argmin}{arg\,min}
+\DeclareMathOperator*{\argmax}{arg\,max}
+
+\DeclarePairedDelimiter\abs{\lvert}{\rvert}
+\DeclarePairedDelimiterX{\inp}[2]{\langle}{\rangle}{#1, #2}
+\DeclarePairedDelimiterX{\infdivx}[2]{(}{)}{%
+  #1\;\delimsize\|\;#2%
+}
+\newcommand{\infdiv}{\infdivx}
+
+\newcommand\iids{\overset{\text{iid}}{\sim}}
+\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
+\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}}
+\newcommand{\Var}{\mathrm{Var}}
+\newcommand{\Cov}{\mathrm{Cov}}
+\newcommand\given[1][]{\:#1\vert\:}
+\newcommand*\mean[1]{\bar{#1}}
+
+% Simple notes for commenting the manuscript
+\usepackage{color}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{document}
+
+\title{Bayesian Continuous Parameters Models}
+
+%\titlerunning{AV uncertainty}
+
+\author{Azat Garifullin}
+
+%
+\maketitle      
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Study diary}
+In his lecture Marko Laine presented a general approach for remote sensing data analysis addressing problems related to big data spatio-temporal analysis with uncertainty quantification. 
+
+Bayesian modelling is a key to solve the problems requiring good priors to model geostatistical data. Gaussian process were presented as a possible prior modelling tool. Due to the nature of the data it is necessary to consider both spatial and temporal dependencies. The lecture considered one of my biggest concerns when dealing with Gaussian processes, namely inversions of huge covariance matrices which is extremely computationally expensive.
+
+Temporal modelling is done using an assumption that spatial and temporal covariance functions are separable what makes the matrices inversion easier due to sparsity induced by the separability. This assumption makes Bayesian filtering and smoothing methods (Kalman filters/smoothers) applicable saving computational resources.
+
+As for spatial modlling Mat\'{e}rn fields were mentioned. Estimation of the parameters of Mat\'{e}rn fields is not an easy task and MCMC is utilised to solve the problem. This also requires the covariance matrices inversions. On one of the first slides Gaussian Markov Fields were mentioned which probably helps to save computational resources.
+
+Nevertheless, even aforementioned methods introducing sparsities seem not to be enough so additional dimensionality reduction methods are employed. The dimensionality reduction methods considered in the lecture are linear what makes it easy to propagate covariance matrices to compressed representations. Kalman filtering and smoothing algorithms can then be trivially extended to work with the compressed representations.
+However, the question of effectiveness of the chosen basis arises. Empirical ortogonal functions and wavelets are a typical choice for the purpose. As a deep learning practioner I cannot to mention that probably self-supervised deep bayesian models could be used for Mat\'{e}rn fields compression. It is well known that neural networks converge to Gaussian processes in a limit of infinite sizes. Therefore, one can probably compress Gaussian processes using Gaussian processes. Gaussian process latent variable models also exist. However, any of the mentioned self-supervised approaches would require to invert the covariance matrices. Thus, it seems that simple linear methods is a best choice at the moment.
+
+Another interesting things are mentioned in the ``Alternative approaches'' slide:
+\begin{enumerate}
+    \item Variational methods with L-BFGS. This made me to remember the talk during the reading group where EnKF inversion was discussed. Nowadays many automatic differentiation frameworks are available. One can probably use that to make efficient MAP estimators. Or to even create new kind of amortised variational approximations, e.g. like Stein Variational Gradient Descent. Still need to invert matrices.
+    \item Transport maps. It reminded me normalising flows, but more scientific.
+\end{enumerate}
+
+The most interesting thing is of course that they actually made things work, whereas we are struggling with getting our chains mixed well for toy problems.
+
+\end{document}
\ No newline at end of file

multiESS @ f64f3102

+Subproject commit f64f3102d404e25c1db704baa30a77d9dc3ada9d

figs/cloudy_matern/mv_ess.txt

+Multivariate ESS: 91.73057089406737
\ No newline at end of file