Commit 92372bb8 authored by LOUIS Maxime's avatar LOUIS Maxime

Adding deep pga. Minor reading/writing improvements

parent d78c70c5
......@@ -69,8 +69,8 @@ class GradientAscent(AbstractEstimator):
# Uncomment for a check of the gradient for the model !
# WARNING: don't forget to comment the update_fixed_effects method of the model !
# print("Checking the model gradient:")
# self._check_model_gradient()
#print("Checking the model gradient:")
#self._check_model_gradient()
self.current_attachment, self.current_regularity, gradient = self._evaluate_model_fit(self.current_parameters,
with_grad=True)
......@@ -262,7 +262,7 @@ class GradientAscent(AbstractEstimator):
attachment, regularity, gradient = self._evaluate_model_fit(self.current_parameters, with_grad=True)
parameters = copy.deepcopy(self.current_parameters)
epsilon = 1e-8
epsilon = 1e-3
for key in gradient.keys():
if key in ['template_data', 'momenta', 'modulation_matrix', 'sources']: continue
......
......@@ -128,8 +128,7 @@ class McmcSaem(AbstractEstimator):
self._maximize_over_fixed_effects()
fixed_effects_after_maximization = self.statistical_model.get_fixed_effects()
fixed_effects = {key: value + step * (fixed_effects_after_maximization[key] - value)
for key, value in fixed_effects_before_maximization.items()}
# fixed_effects = fixed_effects_after_maximization
for key, value in fixed_effects_before_maximization.items()}
self.statistical_model.set_fixed_effects(fixed_effects)
# Averages the random effect realizations in the concentration phase.
......@@ -256,10 +255,12 @@ class McmcSaem(AbstractEstimator):
return aux ** - 0.9
def _initialize_number_of_burn_in_iterations(self):
if self.max_iterations > 4000:
self.number_of_burn_in_iterations = self.max_iterations - 2000
else:
self.number_of_burn_in_iterations = int(self.max_iterations / 2)
if self.number_of_burn_in_iterations is None:
# Because some models will set it manually (e.g. deep Riemannian models)
if self.max_iterations > 4000:
self.number_of_burn_in_iterations = self.max_iterations - 2000
else:
self.number_of_burn_in_iterations = int(self.max_iterations / 2)
def _initialize_acceptance_rate_information(self):
# Initialize average_acceptance_rates.
......
......@@ -80,8 +80,8 @@ class GenericGeodesic:
else:
dt = (self.t0 - self.tmin) / (self.backward_exponential.number_of_time_points - 1)
j = int((time_np - self.tmin) / dt) + 1
assert times[j - 1] <= time_np, "{} {} {}".format(j, time_np, times[j-1])
assert times[j] >= time_np
assert times[j - 1] <= time_np + 1e-3, "{} {} {}".format(j, time_np, times[j-1])
assert times[j] + 1e-3 >= time_np, "{} {} {}".format(j, time_np, times[j])
else:
if self.forward_exponential.number_of_time_points <= 2:
j = len(times) - 1
......@@ -91,7 +91,7 @@ class GenericGeodesic:
int((time_np - self.t0) / dt) + self.backward_exponential.number_of_time_points)
assert times[j - 1] <= time_np
assert times[j] >= time_np
assert times[j] >= time_np, '{}, {}, {}, {}'.format(times[j], time_np, len(times), j)
weight_left = (times[j] - t) / (times[j] - times[j - 1])
weight_right = (t - times[j - 1]) / (times[j] - times[j - 1])
......
......@@ -118,16 +118,16 @@ class ScalarNet(AbstractNet):
super(ScalarNet, self).__init__()
self.layers = nn.ModuleList([nn.Linear(in_dimension, in_dimension),
nn.Tanh(),
nn.Linear(in_dimension, out_dimension, bias=True),
#nn.Linear(in_dimension, out_dimension, bias=True),
#nn.Tanh(),
#nn.Linear(out_dimension, out_dimension, bias=True),
#nn.Tanh(),
nn.Linear(in_dimension, out_dimension),
nn.Tanh(),
nn.Linear(out_dimension, out_dimension, bias=True),
nn.Tanh(),
nn.Linear(out_dimension, out_dimension, bias=True),
nn.Tanh(),
nn.Linear(out_dimension, out_dimension, bias=True),
nn.ELU(),
nn.Linear(out_dimension, out_dimension, bias=True),
nn.ELU()
#nn.Linear(out_dimension, out_dimension),
#nn.Tanh(),
nn.Linear(out_dimension, out_dimension),
nn.Tanh()
])
self.update()
......@@ -180,6 +180,57 @@ class ImageNet2d(AbstractNet):
else:
return x.squeeze(1).squeeze(0)
# This net automatically outputs 128 x 128 images. (to be improved)
class ImageNet2d128(AbstractNet):
def __init__(self, in_dimension=2):
super(ImageNet2d128, self).__init__()
ngf = 2
self.layers = nn.ModuleList([
nn.Linear(in_dimension, in_dimension),
nn.Tanh(), # was added 29/03 12h
nn.ConvTranspose2d(in_dimension, 32 * ngf, 4, stride=4),
# nn.ELU(), # was removed, same day 14h
# nn.ConvTranspose2d(32 * ngf, 16 * ngf, 2, stride=2, bias=False),
nn.ELU(),
# nn.BatchNorm2d(num_features=16*ngf),
nn.ConvTranspose2d(32 * ngf, 16 * ngf, 2, stride=2),
nn.ELU(),
# nn.BatchNorm2d(num_features=8*ngf),
nn.ConvTranspose2d(16 * ngf, 8 * ngf, 2, stride=2),
nn.ELU(),
# nn.BatchNorm2d(num_features=4*ngf),
nn.ConvTranspose2d(8 * ngf, 4 * ngf, 2, stride=2),
nn.ELU(),
# nn.BatchNorm2d(num_features=2*ngf),
nn.ConvTranspose2d(4 * ngf, 2 * ngf, 2, stride=2),
nn.ELU(),
# nn.BatchNorm2d(num_features=2*ngf),
nn.ConvTranspose2d(2 * ngf, 1, 2, stride=2),
nn.ELU()
])
self.update()
def forward(self, x):
a = x.size()
# if len(a) == 2:
# x = x.unsqueeze(2).unsqueeze(3)
# else:
# x = x.unsqueeze(0).unsqueeze(2).unsqueeze(3)
for i, layer in enumerate(self.layers):
if i == 0:
x = layer(x)
if len(a) == 2:
x = x.unsqueeze(2).unsqueeze(3)
else:
x = x.unsqueeze(0).unsqueeze(2).unsqueeze(3)
else:
x = layer(x)
if len(a) == 2:
return x.squeeze(1)
else:
return x.squeeze(1).squeeze(0)
# This net automatically outputs 64 x 64 x 64 images. (to be improved)
class ImageNet3d(AbstractNet):
......@@ -215,3 +266,46 @@ class ImageNet3d(AbstractNet):
else:
return x.squeeze(1).squeeze(0)
# This net automatically outputs 28 x 28 images. (to be improved)
class MnistNet(AbstractNet):
def __init__(self, in_dimension=2):
super(MnistNet, self).__init__()
self.layers = nn.ModuleList([
nn.Linear(in_dimension, 8*2*2),
nn.Tanh(),
nn.ConvTranspose2d(8, 16, 3, stride=2), # b, 16, 5, 5
# nn.ELU(), # was removed, same day 14h
# nn.ConvTranspose2d(32 * ngf, 16 * ngf, 2, stride=2, bias=False),
nn.ELU(),
# nn.BatchNorm2d(num_features=16*ngf),
nn.ConvTranspose2d(16, 8, 5, stride=3, padding=1), # b, 8, 15, 15
nn.ELU(),
# nn.BatchNorm2d(num_features=8*ngf),
nn.ConvTranspose2d(8, 1, 2, stride=2, padding=1), # b, 1, 28, 28
nn.ELU(),
])
self.update()
def forward(self, x):
a = x.size()
# if len(a) == 2:
# x = x.unsqueeze(2).unsqueeze(3)
# else:
# x = x.unsqueeze(0).unsqueeze(2).unsqueeze(3)
for i, layer in enumerate(self.layers):
#print(x.size())
if i == 0:
x = layer(x)
if len(a) == 2:
x = x.unsqueeze(2).unsqueeze(3)
else:
x = x.unsqueeze(0).unsqueeze(2).unsqueeze(3)
x = x.view(len(x), 8, 2, 2)
else:
x = layer(x)
if len(a) == 2:
return x.squeeze(1)
else:
return x.squeeze(1).squeeze(0)
\ No newline at end of file
This diff is collapsed.
......@@ -20,9 +20,10 @@ from pydeformetrica.src.support.probability_distributions.multi_scalar_inverse_w
MultiScalarInverseWishartDistribution
from pydeformetrica.src.support.probability_distributions.multi_scalar_normal_distribution import \
MultiScalarNormalDistribution
from pydeformetrica.src.core.model_tools.manifolds.metric_learning_nets import ScalarNet, ImageNet2d, ImageNet3d
from pydeformetrica.src.core.model_tools.manifolds.metric_learning_nets import ScalarNet, ImageNet2d, ImageNet3d, ImageNet2d128
import matplotlib.pyplot as plt
plt.switch_backend('agg')
from joblib import Parallel, delayed
from torch import nn
from torch import optim
......@@ -313,7 +314,7 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
print("tmin", self.spatiotemporal_reference_frame.geodesic.tmin, "tmax", self.spatiotemporal_reference_frame.geodesic.tmax)
plt.savefig(os.path.join(Settings().output_dir, "Latent_space_coordinates.pdf"))
# plt.savefig(os.path.join(Settings().output_dir, "Latent_space_coordinates.pdf"))
plt.clf()
lsd_observations = torch.from_numpy(lsd_observations).type(Settings().tensor_scalar_type)
......@@ -377,17 +378,18 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
def _individual_RER_to_torch_tensors(self, individual_RER, with_grad):
onset_ages = Variable(torch.from_numpy(individual_RER['onset_age']).type(Settings().tensor_scalar_type),
requires_grad=with_grad)
log_accelerations = Variable(torch.from_numpy(individual_RER['log_acceleration']).type(Settings().tensor_scalar_type),
log_accelerations = Variable(torch.from_numpy(individual_RER['log_acceleration'])
.type(Settings().tensor_scalar_type),
requires_grad=with_grad)
sources = None
if not self.no_parallel_transport:
sources = Variable(torch.from_numpy(individual_RER['sources']),
requires_grad=with_grad)\
.type(Settings().tensor_scalar_type)
sources = Variable(torch.from_numpy(individual_RER['sources']).type(Settings().tensor_scalar_type),
requires_grad=with_grad)
return onset_ages, log_accelerations, sources
return onset_ages, log_accelerations, sources
return onset_ages, log_accelerations, None
def _compute_residuals(self, dataset, v0, p0, metric_parameters, modulation_matrix,
log_accelerations, onset_ages, sources, with_grad=True):
......@@ -401,7 +403,6 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
number_of_subjects = dataset.number_of_subjects
t_begin = time.time()
residuals = []
# # Multi-threaded version
......@@ -414,8 +415,6 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
# args.append((i, j, Settings().serialize(),
# self.spatiotemporal_reference_frame.get_position_exponential(t, sources[i]),
# target))
# t_end_copy = time.time()
# # print("time copying everything", round(1000 * (t_end_copy - t_begin)), "ms")
#
# results = Parallel(n_jobs=Settings().number_of_threads)(delayed(compute_exponential_and_attachment)(arg)
# for arg in args)
......@@ -444,6 +443,7 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
predicted_values_i[j] = self.spatiotemporal_reference_frame.get_position(t)
else:
reference_time = self.get_reference_time()
latent_coordinates_i = Variable(torch.zeros(len(predicted_values_i), self.latent_space_dimension)).type(Settings().tensor_scalar_type)
for j, t in enumerate(absolute_times[i]):
if sources is not None:
......@@ -463,9 +463,6 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
assert Settings().dimension == 3
residuals.append(torch.sum(torch.sum(torch.sum(residuals_i.view(targets_torch.size()), 1), 1, 1)))
# t_end = time.time()
# print("Computing the", dataset.total_number_of_observations,"residuals (after update of the reference frame) took", round(1000*(t_end - t_begin)), "ms")
return residuals
def _compute_individual_attachments(self, residuals):
......@@ -527,12 +524,10 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
if sources is not None:
for i in range(number_of_subjects):
# regularity += torch.norm(sources[i], 1)
regularity += self.individual_random_effects['sources'].compute_log_likelihood_torch(sources[i])
# Noise random effect
regularity -= 0.5 * number_of_subjects \
* math.log(self.fixed_effects['noise_variance'])
regularity -= 0.5 * number_of_subjects * math.log(self.fixed_effects['noise_variance'])
return regularity
......@@ -791,7 +786,7 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
# Set the modulation_matrix prior standard deviation to the deformation kernel width.
self.priors['modulation_matrix'].set_variance_sqrt(1.)
def initialize_deep_metric_learning(self):
def initialize_deep_metric_learning(self, obs_shape=(64,64)):
"""
initialize the neural network and the metric parameters in the model fixed effects
"""
......@@ -802,8 +797,13 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
self.net = ScalarNet(in_dimension=self.latent_space_dimension, out_dimension=Settings().dimension)
elif self.observation_type == 'image':
if Settings().dimension == 2:
print("Defaulting Image net output dimension to 64 x 64")
self.net = ImageNet2d(in_dimension=self.latent_space_dimension)
a, b = obs_shape
if a == 64:
print("Defaulting Image net output dimension to 64 x 64")
self.net = ImageNet2d(in_dimension=self.latent_space_dimension)
elif a == 128:
print("Defaulting Image net output dimension to 64 x 64")
self.net = ImageNet2d128(in_dimension=self.latent_space_dimension)
elif Settings().dimension == 3:
print("Defaulting Image net output dimension to 64 x 64 x 64")
self.net = ImageNet3d(in_dimension=self.latent_space_dimension)
......@@ -874,7 +874,11 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
self._clean_and_create_directory(os.path.join(Settings().output_dir, "geodesic_trajectory"))
self._write_image_trajectory(range(len(times_geodesic)), geodesic_values, os.path.join(Settings().output_dir, "geodesic_trajectory"), "geodesic")
times_parallel_curves = np.linspace(np.min(times_geodesic)+1e-5, np.max(times_geodesic)-1e-5, 20)
if self.observation_type == 'image':
times_parallel_curves = np.linspace(np.min(times_geodesic)+1e-5, np.max(times_geodesic)-1e-5, 20)
else:
times_parallel_curves = np.linspace(np.min(times_geodesic) + 1e-5, np.max(times_geodesic) - 1e-5, 200)
# Saving a txt file with the trajectory.
write_2D_array(times_geodesic, self.name + "_reference_geodesic_trajectory_times.txt")
......@@ -915,6 +919,7 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
if self.observation_type == 'scalar':
write_2D_array(trajectory_pos, self.name+"_source_" + str(i) + "_pos.txt")
write_2D_array(trajectory_neg, self.name+"_source_" + str(i) + "_neg.txt")
write_2D_array(times_parallel_curves, self.name + "_times_parallel_curves.txt")
self._plot_scalar_trajectory(times_geodesic, geodesic_values)
self._plot_scalar_trajectory(times_parallel_curves, trajectory_pos, linestyles=['dashed'] * len(trajectory_pos[0]))
self._plot_scalar_trajectory(times_parallel_curves, trajectory_neg, linestyles=['dotted'] * len(trajectory_neg[0]))
......@@ -965,7 +970,7 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
nb_plot_to_make = 10000
if sample:
nb_plot_to_make = float("inf")
subjects_per_plot = 1
subjects_per_plot = 3
predictions = []
subject_ids = []
times = []
......@@ -1042,7 +1047,9 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
else:
self._clean_and_create_directory(os.path.join(Settings().output_dir, "subject_"+str(i)))
if Settings().dimension == 2:
trajectory = trajectory.reshape(len(trajectory), 64, 64)
if len(trajectory.shape) < 3:
image_dimension = int(math.sqrt(len(trajectory[0])))
trajectory = trajectory.reshape(len(trajectory), image_dimension, image_dimension)
elif Settings().dimension == 3:
trajectory = trajectory.reshape(len(trajectory), 64, 64, 64)
self._write_image_trajectory(dataset.times[i], trajectory,
......@@ -1126,7 +1133,7 @@ class LongitudinalMetricLearning(AbstractStatisticalModel):
lsd_observations = self._get_lsd_observations(individual_RER, dataset)
write_2D_array(lsd_observations, self.name + '_latent_space_positions.txt')
plt.scatter(lsd_observations[:, 0], lsd_observations[:, 1])
plt.savefig(os.path.join(Settings().output_dir, self.name+'_all_latent_space_positions.pdf'))
# plt.savefig(os.path.join(Settings().output_dir, self.name+'_all_latent_space_positions.pdf'))
plt.clf()
import numpy as np
class LongitudinalDataset:
"""
......@@ -54,7 +57,21 @@ class LongitudinalDataset:
shape = img.get_points().shape
else:
assert img.get_points().shape == shape,\
"Different images dimensions were detected... please code the projection or normalize before"
"Different images dimensions were detected."
def order_observations(self):
""" sort the visits for each individual, by time"""
for i in range(len(self.times)):
arg_sorted_times = np.argsort(self.times[i])
sorted_times = np.sort(self.times[i])
sorted_deformable_objects = []
for j in arg_sorted_times:
sorted_deformable_objects.append(self.deformable_objects[i][j])
self.times[i] = sorted_times
self.deformable_objects[i] = sorted_deformable_objects
......
......@@ -5,7 +5,6 @@ import numpy as np
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + os.path.sep + '../../../')
from pydeformetrica.src.support.utilities.general_settings import Settings
from vtk import vtkPolyDataWriter, vtkPolyData, vtkPoints, vtkDoubleArray
from scipy.misc import toimage, imsave
import nibabel as nib
......@@ -134,20 +133,27 @@ def write_control_points_and_momenta_vtk(control_points, momenta, name):
writer.Update()
def write_2d_image(img_data, name):
def write_2d_image(img_data, name, fmt='.npy'):
"""
img_data is a (pixels * pixels) np array
"""
# Watch out ! only images with npy format (numpy) are used in here !
# imsave(os.path.join(Settings().output_dir, name), img_data)
img = toimage(img_data)
if name.find(Settings().output_dir+"/") >= 0:
# img.save(name)
np.save(name, img_data)
if fmt == '.npy':
img = toimage(img_data)
if name.find(Settings().output_dir+"/") >= 0:
# img.save(name)
np.save(name, img_data)
else:
# img.save(os.path.join(Settings().output_dir, name))
np.save(os.path.join(Settings().output_dir, name), img_data)
else:
# img.save(os.path.join(Settings().output_dir, name))
np.save(os.path.join(Settings().output_dir, name), img_data)
img = toimage(img_data)
if name.find(Settings().output_dir + "/") >= 0:
img.save(name)
else:
img.save(os.path.join(Settings().output_dir, name))
def write_3d_image(img_data, name):
......
......@@ -65,7 +65,7 @@ def create_scalar_dataset(group, observations, timepoints):
for i in range(len(observations)):
if group[i] == subject_id:
times_subject.append(timepoints[i])
scalars_subject.append(observations[i].get_points())
scalars_subject.append(observations[i])
assert len(times_subject) > 0, subject_id
assert len(times_subject) == len(scalars_subject)
times.append(np.array(times_subject))
......@@ -105,6 +105,8 @@ def create_image_dataset(group, observations, timepoints):
longitudinal_dataset.deformable_objects = images
longitudinal_dataset.number_of_subjects = len(subject_ids)
longitudinal_dataset.total_number_of_observations = len(timepoints)
longitudinal_dataset.check_image_shapes()
longitudinal_dataset.order_observations()
return longitudinal_dataset
......@@ -147,6 +149,7 @@ def read_and_create_image_dataset(dataset_filenames, visit_ages, subject_ids, te
longitudinal_dataset.deformable_objects = deformable_objects_dataset
longitudinal_dataset.update()
longitudinal_dataset.check_image_shapes()
longitudinal_dataset.order_observations()
return longitudinal_dataset
......
......@@ -48,13 +48,14 @@ class DeformableObjectReader:
out_object.set_connectivity(connectivity)
elif object_type.lower() == 'PointCloud'.lower():
out_object = PointCloud()
points = DeformableObjectReader.read_vtk_file(object_filename, extract_connectivity=False)
out_object.set_points(points)
elif object_type.lower() == 'Landmark'.lower():
out_object = Landmark()
points = DeformableObjectReader.read_vtk_file(object_filename, extract_connectivity=False)
out_object.set_points(self._extract_points(poly_data))
out_object.set_points(points)
out_object.update()
......@@ -146,25 +147,3 @@ class DeformableObjectReader:
else:
return points
# @staticmethod
# def _extract_points(poly_data):
# number_of_points = poly_data.GetNumberOfPoints()
# dimension = Settings().dimension
# points = np.zeros((number_of_points, dimension))
# for k in range(number_of_points):
# p = poly_data.GetPoint(k)
# points[k, :] = p[0:dimension]
# return points
#
# @staticmethod
# def _extract_connectivity(poly_data, nb_points_per_face):
# """
# extract the list of faces from a poly data, and returns it as a numpy array
# nb_points_per_face is the number of points on each face
# """
# connectivity = np.zeros((poly_data.GetNumberOfCells(), nb_points_per_face))
# for i in range(poly_data.GetNumberOfCells()):
# for j in range(nb_points_per_face):
# connectivity[i, j] = poly_data.GetCell(i).GetPointId(j)
# return connectivity
import os.path
import sys
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + os.path.sep + '../../../')
import torch
import torchvision
from torch.autograd import Variable
from torch import nn
from torch import optim
from torch.utils.data import TensorDataset, DataLoader
import numpy as np
import warnings
import time
# Estimators
from pydeformetrica.src.core.estimators.mcmc_saem import McmcSaem
from pydeformetrica.src.core.estimators.scipy_optimize import ScipyOptimize
from pydeformetrica.src.core.estimators.gradient_ascent import GradientAscent
from pydeformetrica.src.core.estimator_tools.samplers.srw_mhwg_sampler import SrwMhwgSampler
from pydeformetrica.src.support.utilities.general_settings import Settings
from pydeformetrica.src.support.probability_distributions.multi_scalar_normal_distribution import MultiScalarNormalDistribution
from pydeformetrica.src.core.observations.datasets.longitudinal_dataset import LongitudinalDataset
from pydeformetrica.src.core.observations.manifold_observations.image import Image
from pydeformetrica.src.core.model_tools.manifolds.metric_learning_nets import MnistNet
from pydeformetrica.src.core.models.deep_pga import DeepPga
from sklearn.metrics import mean_squared_error
from sklearn.decomposition import PCA
# Runs the deep pga model on 1000 randomly extracted digits for mnist, for varying latent space dimensions.
# Experimental for now.
#for latent_space_dimension in range(2, 20):
for latent_space_dimension in [2]:
Settings().dimension = 2
Settings().output_dir = 'output_' + str(latent_space_dimension)
if not os.path.exists(Settings().output_dir): os.makedirs(Settings().output_dir)
trainset = torchvision.datasets.MNIST('/tmp', train=True, download=True)
train_dataset = LongitudinalDataset()
test_dataset = LongitudinalDataset()
train_labels = []
test_labels = []
for elt in trainset:
image, label = elt
#if label == 2:
if len(train_dataset.deformable_objects) < 800:
img = Image()
img.set_points(np.array(image, dtype=float)/np.max(image))
img.update()
train_dataset.deformable_objects.append([img])
train_dataset.subject_ids.append(len(train_dataset.deformable_objects))
train_labels.append(label)
elif len(test_dataset.deformable_objects) < 200:
img = Image()
img.set_points(np.array(image, dtype=float)/np.max(image))
img.update()
test_dataset.deformable_objects.append([img])
test_dataset.subject_ids.append(len(test_dataset.deformable_objects))
test_labels.append(label)
np.savetxt(os.path.join(Settings().output_dir, 'labels_train.txt'), np.array(train_labels))
np.savetxt(os.path.join(Settings().output_dir, 'labels_test.txt'), np.array(test_labels))
a, b = np.array(image).shape
train_dataset.update()
test_dataset.update()
images_data = np.array([elt[0].get_points() for elt in train_dataset.deformable_objects])
images_data_torch = torch.from_numpy(images_data).type(Settings().tensor_scalar_type)
pca = PCA(n_components=latent_space_dimension)
latent_space_positions = pca.fit_transform([elt[0].get_points().flatten() for elt in train_dataset.deformable_objects])
reconstructed = pca.inverse_transform(latent_space_positions)
latent_test = pca.transform([elt[0].get_points().flatten() for elt in test_dataset.deformable_objects])
reconstructed_test = pca.inverse_transform(latent_test)
reconstruction_error_train = mean_squared_error(reconstructed, [elt[0].get_points().flatten() for elt in train_dataset.deformable_objects])
reconstruction_error_test = mean_squared_error(reconstructed_test, [elt[0].get_points().flatten() for elt in test_dataset.deformable_objects])
print('PCA mse on train:', reconstruction_error_train)
print('PCA mse on test:', reconstruction_error_train)
to_write = np.array([reconstruction_error_train, reconstruction_error_test])
np.savetxt(os.path.join(Settings().output_dir, 'pca_reconstruction_error.txt'), to_write)
# We now normalize every latent_space_positions
for i in range(latent_space_dimension):
latent_space_positions[:, i] = (latent_space_positions[:, i] - np.mean(latent_space_positions[:, i]))/np.std(latent_space_positions[:, i])
latent_space_positions_torch = torch.from_numpy(latent_space_positions).type(Settings().tensor_scalar_type)
# We now instantiate the neural network
net = MnistNet(in_dimension=latent_space_dimension)
#net.double()
criterion = nn.MSELoss()
nb_epochs = 50
optimizer = optim.Adam(net.parameters(), lr=1e-3, weight_decay=0)
train_dataset_fit = TensorDataset(latent_space_positions_torch, images_data_torch)
train_dataloader = DataLoader(train_dataset_fit, batch_size=10, shuffle=True)
for epoch in range(nb_epochs):
train_loss = 0
nb_train_batches = 0
for (z, y) in train_dataloader:
nb_train_batches += 1
var_z = Variable(z)
var_y = Variable(y)
predicted = net(var_z)
loss = criterion(predicted, var_y)
net.zero_grad()
loss.backward()
train_loss += loss.cpu().data.numpy()[0]
optimizer.step()
train_loss /= nb_train_batches
print("Epoch {}/{}".format(epoch, nb_epochs),
"Train loss:", train_loss)
metric_parameters = net.get_parameters()
# We instantiate the model
model = DeepPga()
model.net = net
model.set_metric_parameters(metric_parameters)
model.latent_space_dimension = latent_space_dimension
model.set_noise_variance(train_loss / (a*b))
model.update()
model.individual_random_effects['latent_position'].mean = np.zeros((latent_space_dimension,))
model.individual_random_effects['latent_position'].set_variance(1.0)
individual_RER = {}
individual_RER['latent_position'] = latent_space_positions
sampler = SrwMhwgSampler()
estimator = McmcSaem()
estimator.sampler = sampler
# latent positions proposal:
latent_position_proposal_distribution = MultiScalarNormalDistribution()
latent_position_proposal_distribution.set_variance_sqrt(0.1)
sampler.individual_proposal_distributions['latent_position'] = latent_position_proposal_distribution
estimator.sample_every_n_mcmc_iters = 10
estimator.max_iterations = 200
estimator.number_of_burn_in_iterations = 200
estimator.max_line_search_iterations = 10
estimator.convergence_tolerance = 1e-3
estimator.print_every_n_iters = 1
estimator.save_every_n_iters = 30
estimator.dataset = train_dataset
estimator.statistical_model = model
# Initial random effects realizations
estimator.individual_RER = individual_RER
if not os.path.exists(Settings().output_dir): os.makedirs(Settings().output_dir)
model.name = 'DeepPgaModel'
print('')
print('[ update method of the ' + estimator.name + ' optimizer ]')