Commit c817d22f authored by Azat Garifullin's avatar Azat Garifullin
Browse files

gibbs sampling

parents
# Default ignored files
/shelf/
/workspace.xml
gp_mix
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$/../gp_mix" />
<orderEntry type="jdk" jdkName="Python 3.8 (calib) (2)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
</module>
\ No newline at end of file
<component name="InspectionProjectProfileManager">
<profile version="1.0">
<option name="myName" value="Project Default" />
<inspection_tool class="PyPep8NamingInspection" enabled="true" level="WEAK WARNING" enabled_by_default="true">
<option name="ignoredErrors">
<list>
<option value="N803" />
<option value="N806" />
</list>
</option>
</inspection_tool>
</profile>
</component>
\ No newline at end of file
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (calib) (2)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/../gp_mix/.idea/gp_mix.iml" filepath="$PROJECT_DIR$/../gp_mix/.idea/gp_mix.iml" />
</modules>
</component>
</project>
\ No newline at end of file
.idea
__pycache__/
import itertools
import numpy as np
from scipy import spatial
from scipy.special import kv, gamma
def log_mvnpdf(x, mu, sigma):
cov = (sigma**2) * np.eye(mu.size) if np.isscalar(sigma) else (sigma**2)
xi = x - mu
md = 0.5 * np.dot(xi, np.dot(np.linalg.pinv(cov), xi))
sd, ld = np.linalg.slogdet(cov)
nc = 0.5 * (mu.size * np.log(2 * np.pi)) + sd * ld
return -md - nc
def matern_cov(xs, ys, nu, sigma, l):
return matern_cov_from_distances(dist_matrix(xs, ys), nu, sigma, l)
def matern_cov_from_distances(dist_mtx, nu, sigma, l):
s = (sigma**2) * np.power(2.0, 1 - nu) / gamma(nu)
safe_ds = np.maximum(dist_mtx, np.sqrt(np.finfo(np.float32).eps))
ld = np.sqrt(2 * nu) * safe_ds / l
kd = kv(nu, ld)
cov_mtx = s * np.power(ld, nu) * kd
return cov_mtx
def dist_matrix(xs, ys):
xv, yv = np.meshgrid(xs, ys)
xy_grid = np.stack([xv.ravel(), yv.ravel()], axis=1)
dist_mtx = spatial.distance.cdist(xy_grid, xy_grid)
return dist_mtx
def sample_matern_cond_posterior(mean_field, cov_mtx, noise_sigma, fix_mean=True):
i_mtx = np.eye(cov_mtx.shape[0])
nc_inv = np.linalg.inv(cov_mtx + (noise_sigma**2) * i_mtx)
cnc = np.matmul(cov_mtx, nc_inv)
mu = np.matmul(cnc, mean_field)
delta = 0
if not fix_mean:
cov = cov_mtx - np.matmul(cnc, cov_mtx)
chol_cov = np.linalg.cholesky(cov)
n = np.random.normal(size=(cov_mtx.shape[0],))
delta = np.matmul(chol_cov, n)
return mu + delta
def matern_parameter_nlpost(observed, field, field_cov, noise_sigma):
inv_noise_cov = np.eye(np.size(field)) * (1.0 / noise_sigma**2)
field_diff = observed - field
nll = np.matmul(np.matmul(field_diff, inv_noise_cov), field_diff)
field_inv_cov = np.linalg.inv(field_cov)
field_nl_prior = np.matmul(np.matmul(field.T, field_inv_cov), field)
sgn, val = np.linalg.slogdet(field_cov)
log_det_cov = sgn * val
return nll + field_nl_prior + log_det_cov
def matern_nl_marginal_likelihood(observed_field, nu, sigma, l, dist_mtx, noise_sigma):
cov_mtx = matern_cov_from_distances(dist_mtx, nu, sigma, l)
i_mtx = np.eye(cov_mtx.shape[0])
cov = cov_mtx + (noise_sigma**2) * i_mtx
inv_cov = np.linalg.inv(cov)
sgn, val = np.linalg.slogdet(cov)
log_det_cov = sgn * val
return log_det_cov + 0.5 * np.matmul(np.matmul(observed_field, inv_cov), observed_field)
def sample_matern_field(xs, ys, nu, sigma, l):
cov_mtx = matern_cov(xs, ys, nu, sigma, l)
chol_cov = np.linalg.cholesky(cov_mtx)
std_noise = np.random.normal(size=(chol_cov.shape[0],))
return np.matmul(chol_cov, std_noise)
def gibbs_within_metropolis(initial_theta, observed_field, noise_sigma,
proposal_cov, lower, upper,
dist_mtx, initial_field=None,
adaptation_start=None,
n_samples=10000, n_save=1000):
chain = []
nlls = []
p_accs = []
state = np.log(initial_theta)
field = initial_field
if field is None:
field = np.copy(observed_field)
nu, sigma, l = initial_theta
field_cov = matern_cov_from_distances(dist_mtx, nu, sigma, l)
state_nlp = matern_parameter_nlpost(observed_field, field,
field_cov, noise_sigma)
p_mean_intext = np.copy(state)
sd = 5.76 / state.size
eps = np.finfo(np.float32).eps
eps_eye = eps * np.eye(state.size)
proposal_cov += eps_eye
if adaptation_start is None:
adaptation_start = n_samples
for step in range(n_samples):
pr_chol = np.linalg.cholesky(proposal_cov)
dstate = np.matmul(pr_chol, np.random.normal(size=state.shape))
proposal_state = state + dstate
proposal_theta = np.exp(proposal_state)
up_lower = np.all(proposal_theta >= lower)
lo_upper = np.all(proposal_theta <= upper)
if not (lo_upper and up_lower):
print(
"Gibbs-Metropolis [{}]: rejected with out of bounds".format(step))
continue
nu, sigma, l = proposal_theta
proposal_field_cov = matern_cov_from_distances(dist_mtx, nu, sigma, l)
proposal_field = sample_matern_cond_posterior(
observed_field, proposal_field_cov, noise_sigma)
proposal_nlp = matern_parameter_nlpost(
observed_field, proposal_field, proposal_field_cov, noise_sigma)
acceptance = np.exp(state_nlp - proposal_nlp)
accepted = np.random.random() < acceptance
if accepted:
state = proposal_state
field = proposal_field
state_nlp = proposal_nlp
print(
"Gibbs-Metropolis [{}]: acceptance with p={}".format(step, acceptance))
else:
print(
"Gibbs-Metropolis [{}]: rejected with p={}".format(step, acceptance))
if step > (n_samples - n_save):
chain.append((np.exp(state), field))
n_mean_intext = p_mean_intext + (state - p_mean_intext) / (step + 1)
if step >= adaptation_start:
base = (step - 1) / step * proposal_cov
pmmt = np.outer(p_mean_intext, p_mean_intext) * step
nmmt = np.outer(n_mean_intext, n_mean_intext) * (step + 1)
sst = np.outer(state, state)
upd_cov = pmmt - nmmt + sst + eps_eye
proposal_cov = base + sd / step * upd_cov
p_accs.append((acceptance, accepted))
nlls.append(state_nlp)
return chain, p_accs, nlls
import itertools
import numpy as np
from scipy import spatial
from scipy import stats
from scipy import special
from gp_mixes import gp
class MaternParams:
def __init__(self):
pass
def unpack_params(theta):
nu_1, sigma_1, l_1 = theta[:3]
nu_2, sigma_2, l_2 = theta[3:6]
nn_ws = theta[6:]
return nu_1, sigma_1, l_1, nu_2, sigma_2, l_2, nn_ws
def pack_params(nu_1, sigma_1, l_1,
nu_2, sigma_2, l_2,
nn_ws=None):
ws = nn_ws.tolist() if nn_ws is not None else []
return np.array([nu_1, sigma_1, l_1, nu_2, sigma_2, l_2] + ws)
def gating_loglikelihood(field, dist_mtx, nu, sigma, l, noise_Sigma=0):
full_cov = gp.matern_cov_from_distances(dist_mtx, nu, sigma, l) + \
(noise_Sigma ** 2) * np.eye(field.size)
log_likelihood = np.zeros_like(field)
for i in range(field.size):
ks = [k for k in range(field.size) if k != i]
cov_kj = full_cov[np.ix_(ks, ks)]
cov_ij = full_cov[i, ks]
cov_ii = full_cov[i, i]
field_wi = field[ks]
cc = np.matmul(cov_ij.T, cov_kj)
var = np.maximum(cov_ii - np.matmul(cc, cov_ij), np.finfo(np.float32).eps)
mu = np.matmul(cc, field_wi)
log_likelihood[i] = gp.log_mvnpdf(0, mu, np.sqrt(var))
return log_likelihood - special.logsumexp(log_likelihood)
def gating_loglikelihood_patch(field, dist_mtx, nu, sigma, l,
n_x, n_y, patch_size=2):
full_cov = gp.matern_cov_from_distances(dist_mtx, nu, sigma, l)
log_likelihood = np.zeros_like(field)
indices = itertools.product(range(0, n_y, patch_size), range(0, n_x, patch_size))
for i, j in indices:
patch_grid = itertools.product(range(i, i + patch_size), range(j, j + patch_size))
patch_indices = [ip * n_x + jp for ip, jp in patch_grid]
rest_indices = [k for k in range(field.size) if not (k in patch_indices)]
cov_kj = full_cov[np.ix_(rest_indices, rest_indices)]
cov_ij = full_cov[np.ix_(rest_indices, patch_indices)]
cov_ii = full_cov[np.ix_(patch_indices, patch_indices)]
field_wi = field[rest_indices]
cc = np.matmul(cov_ij.T, cov_kj)
cov_patch = cov_ii - np.matmul(cc, cov_ij)
mu = np.matmul(cc, field_wi)
log_likelihood[patch_indices] = gp.log_mvnpdf(0, mu, cov_patch)
return log_likelihood - special.logsumexp(log_likelihood)
import pickle
import random
import numpy as np
from scipy import optimize
from scipy import spatial
from sklearn.decomposition import PCA
from sklearn.mixture import GaussianMixture
from sklearn.neural_network import MLPClassifier
import matplotlib
from matplotlib import pyplot as plt
from gp_mixes import gp
from gp_mixes import utils
from gp_mixes import mix_field as mf
from gp_mixes import utils
def load_ground_truth(path):
with open(path, 'rb') as gt_file:
simulation = pickle.load(gt_file)
return simulation
def make_gp_loss(field, dist_mtx, noise_sigma):
def gp_loss(theta):
nu, sigma, l = theta
return gp.matern_nl_marginal_likelihood(field, nu, sigma, l,
dist_mtx, noise_sigma)
return gp_loss
def main():
simulation = load_ground_truth('ground_truth.pickle')
gt_theta, gt_mix_field, mix_field, noise_sigma, xs, ys = simulation
gabor_bank = utils.make_gabor_bank()
responses = utils.apply_bank_of_filters(mix_field.reshape((ys.size, xs.size)),
gabor_bank)
texture_features = responses.reshape(mix_field.shape + (responses.shape[-1],))
clustering = GaussianMixture(n_components=2)
clustering.fit(texture_features)
labels = clustering.predict(texture_features)
theta_bounds = [(0.5, 3.5), (0.5, 3.5), (0.5, 3.5)]
xv, yv = np.meshgrid(xs, ys)
xy_grid = np.stack([xv.ravel(), yv.ravel()], axis=1)
# xy_grid_1 = xy_grid[labels == 0]
# dist_mtx_1 = spatial.distance.cdist(xy_grid_1, xy_grid_1)
# loss_1 = make_gp_loss(mix_field[labels == 0], dist_mtx_1, noise_sigma)
# print("debug: GT loss_1: {}".format(loss_1(gt_theta[:3])))
# print("debug: GT_2 loss_1: {}".format(loss_1(gt_theta[3:])))
# opt_res_1 = optimize.differential_evolution(loss_1, bounds=theta_bounds, tol=1e-4,
# disp=True, seed=42)
# print(opt_res_1.x)
#
# xy_grid_2 = xy_grid[labels == 1]
# dist_mtx_2 = spatial.distance.cdist(xy_grid_2, xy_grid_2)
# loss_2 = make_gp_loss(mix_field[labels == 1], dist_mtx_2, noise_sigma)
# print("debug: GT loss_2: {}".format(loss_2((gt_theta[3:]))))
# print("debug: GT_2 loss_1: {}".format(loss_2((gt_theta[:3]))))
# opt_res_2 = optimize.differential_evolution(loss_2, bounds=theta_bounds, tol=1e-4,
# disp=True, seed=42)
# print(opt_res_2.x)
gating = MLPClassifier(hidden_layer_sizes=(16, 16), solver='lbfgs', max_iter=300)
gating.fit(xy_grid, labels)
gate_labels = gating.predict(xy_grid)
gate_probas = gating.predict_proba(xy_grid)
# with open('optimization.pickle', 'wb') as optfile:
# thetas = (opt_res_1.x, opt_res_2.x)
# pickle.dump((thetas, gating), optfile)
plt.figure()
plt.imshow(labels.reshape(ys.size, xs.size), cmap='gray')
plt.figure()
plt.imshow(gate_labels.reshape(ys.size, xs.size), cmap='gray')
plt.figure()
plt.imshow(gate_probas[:, 1].reshape(ys.size, xs.size), cmap='gray')
pca = PCA(n_components=1)
r_texture_fs = pca.fit_transform(texture_features)
plt.figure()
plt.imshow(r_texture_fs.reshape(ys.size, xs.size))
plt.figure()
plt.imshow(mix_field.reshape(ys.size, xs.size))
dist_mtx = gp.dist_matrix(xs, ys)
gl_1 = mf.gating_loglikelihood_patch(gt_mix_field, dist_mtx, gt_theta[0], gt_theta[1], gt_theta[2],
n_x=xs.size, n_y=ys.size)
gl_2 = mf.gating_loglikelihood_patch(gt_mix_field, dist_mtx, gt_theta[3], gt_theta[4], gt_theta[5],
n_x=xs.size, n_y=ys.size)
gl_cat = np.stack([gl_1, gl_2], axis=1)
mask = np.argmax(gl_cat, axis=1)
plt.figure()
plt.imshow(mask.reshape((24, 24)), cmap='gray')
plt.figure()
plt.imshow(gl_1.reshape((24, 24)), cmap='gray')
plt.figure()
plt.imshow(gl_2.reshape((24, 24)), cmap='gray')
plt.show()
if __name__ == '__main__':
random.seed(42)
np.random.seed(42)
main()
\ No newline at end of file
import pickle
import random
import numpy as np
from scipy import optimize
from scipy import spatial
from scipy import special
from sklearn.decomposition import PCA
from sklearn.mixture import GaussianMixture
from sklearn.neural_network import MLPClassifier
import matplotlib
from matplotlib import pyplot as plt
from gp_mixes import gp
from gp_mixes import utils
from gp_mixes import mix_field as mf
from gp_mixes import utils
def load_ground_truth(path):
with open(path, 'rb') as gt_file:
simulation = pickle.load(gt_file)
return simulation
def load_optimization(path):
with open(path, 'rb') as opt_file:
thetas, gating = pickle.load(opt_file)
return thetas, gating
def main():
gt = load_ground_truth('ground_truth.pickle')
gt_theta, gt_mix_field, mix_field, noise_sigma, xs, ys = gt
thetas, gating = load_optimization('optimization.pickle')
theta_1, theta_2 = thetas
xv, yv = np.meshgrid(xs, ys)
xy_grid = np.stack([xv.ravel(), yv.ravel()], axis=1)
dist_mtx = spatial.distance.cdist(xy_grid, xy_grid)
theta_lower = [0.5, 0.5, 0.5]
theta_upper = [3.5, 3.5, 3.5]
proposal_1 = 1e-3 * np.eye(3, dtype=np.float32)
proposal_2 = 1e-3 * np.eye(3, dtype=np.float32)
total_chain_1 = []
total_chain_2 = []
all_fields = []
all_masks = []
n_iters = 200
for iteration in range(n_iters):
gate_labels = gating.predict(xy_grid)
all_masks.append(gate_labels)
new_field = np.zeros_like(mix_field)
gate_labels_1 = gate_labels == 0
if np.count_nonzero(gate_labels_1) != 0:
xy_grid_1 = xy_grid[gate_labels_1]
dist_mtx_1 = spatial.distance.cdist(xy_grid_1, xy_grid_1)
chain_1, _, _ = gp.gibbs_within_metropolis(
theta_1, mix_field[gate_labels_1], noise_sigma, proposal_1,
theta_lower, theta_upper, dist_mtx_1,
adaptation_start=10,
n_samples=100, n_save=100)
thetas_chain_1, fields_1 = zip(*chain_1)
proposal_1 = np.cov(np.log(thetas_chain_1).T)
total_chain_1 += thetas_chain_1
new_field[gate_labels == 0] = np.mean(fields_1, axis=0)
theta_1 = thetas_chain_1[-1]
else:
print("Skipping label 0")
gate_labels_2 = gate_labels == 1
if np.count_nonzero(gate_labels_2) != 0:
xy_grid_2 = xy_grid[gate_labels_2]
dist_mtx_2 = spatial.distance.cdist(xy_grid_2, xy_grid_2)
chain_2, _, _ = gp.gibbs_within_metropolis(
theta_2, mix_field[gate_labels_2], noise_sigma, proposal_2,
theta_lower, theta_upper, dist_mtx_2,
adaptation_start=10,
n_samples=100, n_save=100)
thetas_chain_2, fields_2 = zip(*chain_2)
proposal_2 = np.cov(np.log(thetas_chain_2).T)
total_chain_2 += thetas_chain_2
new_field[gate_labels == 1] = np.mean(fields_2, axis=0)
theta_2 = thetas_chain_2[-1]
else:
print("Skipping label 1")
all_fields.append(new_field)
gl_1 = mf.gating_loglikelihood(mix_field, dist_mtx,
theta_1[0], theta_1[1], theta_1[2],
noise_sigma)
gl_2 = mf.gating_loglikelihood(mix_field, dist_mtx,
theta_2[0], theta_2[1], theta_2[2],
noise_sigma)
upd_ll = np.stack([gl_1, gl_2], axis=1)
prev_ll = gating.predict_log_proba(xy_grid)
prev_ll -= special.logsumexp(prev_ll)
gating_ll = prev_ll + upd_ll
gating_max = np.argmax(gating_ll, axis=1)
gating.fit(xy_grid, gating_max)
with open('mcmc.pickle', 'wb') as mcmc_file:
mcmc_res = (total_chain_1, total_chain_2,
all_fields, all_masks)
pickle.dump(mcmc_res, mcmc_file)
if __name__ == '__main__':
random.seed(42)
np.random.seed(42)
main()
import random
import pickle
import cv2
import numpy as np
from scipy import optimize
from sklearn.metrics import log_loss
import matplotlib
from matplotlib import pyplot as plt
import gp_mixes.mix_field
from gp_mixes import gp
from gp_mixes import mix_field as mf
def main():
mask = cv2.imread('field_mask.bmp', cv2.IMREAD_GRAYSCALE)
mask = mask >= 127
n_x = 24
xs = np.linspace(0, 10, n_x)
n_y = 24
ys = np.linspace(0, 10, n_y)
fst_nu = 1.0