Commit 1a77bf6c authored by Alexandre Bône's avatar Alexandre Bône
parents 9dc77e33 926455a5
Pipeline #57435700 failed with stages
in 25 minutes and 34 seconds
......@@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
- Improve TorchKernel
- Upgrades: Closes #45
- pytorch to 1.0.1
- pykeops to 1.0
- pykeops to 1.0.1
- numpy to 1.16
- nibabel to 2.3.3
- matplotlib to 3.0
......@@ -34,6 +34,8 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
- Refactoring of the deformable object class, now more robust.
- New VTK reader, based on the conda-available VTK library (v8). Closes #50.
- Generalized and harmonized use of the logger, along with automatic txt file dump. Closes #54.
- Renamed `number-of-threads` to `number-of-processes` to better reflect reality
- Replace `use_cuda` boolean with `gpu_mode` enum. Available gpu_modes are: auto, full, none, kernel
## [4.1.0] - 2018-09-14
### Added
......
......@@ -17,16 +17,24 @@ build:
requirements:
build:
- python
# - numpy # needed because of pykeops 1.0 bug
- setuptools
- pip
run:
- python {{ python }}
- numpy=1.16.2
- scikit-learn=0.20.3
- matplotlib=2.2.2
- nibabel=2.3.3
- pillow=5.4.1
- pytorch=1.0.1
- cudatoolkit=8.0
- psutil=5.4.8
- vtk=8.2.0
# {% if BUILD_GUI == "True" %}
# {% if data.get('build_gui') == "True" %}
- pyqt=5.9
# {% endif %}
test:
commands:
......
......@@ -5,17 +5,20 @@ channels:
- anaconda
dependencies:
- python=3.7
- numpy=1.16.2
- scikit-learn=0.20.3
- matplotlib=2.2.2
- nibabel=2.3.3
- pillow=5.1.0
- h5py=2.9.0 # fix: h5py conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated
- pytorch=1.0.1
- vtk==8.2.0
- cudatoolkit=8.0
- pyqt=5.9
- vtk=8.2.0
- psutil=5.4.8
- ninja=1.9
- pip
- pip:
- PyQt5>=5.11
- gputil==1.3.0
- pykeops==1.0
- pykeops==1.0.1
- cmake>=3.10
- sklearn
- scikit-learn
......@@ -71,14 +71,11 @@ def build_deformetrica():
'Topic :: Software Development :: Libraries'
],
install_requires=[
'numpy>=1.10',
'h5py>=2.8', # fix: h5py conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated
'pykeops==1.0.1',
'gputil>=1.3',
'pykeops==1.0',
'scikit-learn==0.20.3',
'PyQt5>=5.11',
# 'PyQt5>=5.11',
],
# extra_link_args=['-headerpad_max_install_names']
extra_link_args=['-Wl,-headerpad_max_install_names']
)
......@@ -118,13 +115,10 @@ def build_deformetrica_nox():
'Topic :: Software Development :: Libraries'
],
install_requires=[
'numpy>=1.10',
'h5py>=2.8', # fix: h5py conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated
'pykeops==1.0.1',
'gputil>=1.3',
'pykeops==1.0',
'scikit-learn==0.20.3',
],
# extra_link_args=['-headerpad_max_install_names']
extra_link_args=['-Wl,-headerpad_max_install_names']
)
......
......@@ -9,7 +9,7 @@ import time
import torch
import numpy as np
from core import default
from core import default, GpuMode
from core.estimators.gradient_ascent import GradientAscent
from core.estimators.mcmc_saem import McmcSaem
from core.estimators.scipy_optimize import ScipyOptimize
......@@ -553,10 +553,11 @@ class Deformetrica:
# Initializes variables that will be checked.
#
if estimator_options is not None:
if 'use_cuda' not in estimator_options:
estimator_options['use_cuda'] = default.use_cuda
else:
default.update_use_cuda(estimator_options['use_cuda'])
if 'gpu_mode' not in estimator_options:
estimator_options['gpu_mode'] = default.gpu_mode
if estimator_options['gpu_mode'] is GpuMode.FULL and not torch.cuda.is_available():
logger.warning("GPU computation is not available, falling-back to CPU.")
estimator_options['gpu_mode'] = GpuMode.NONE
if 'state_file' not in estimator_options:
estimator_options['state_file'] = default.state_file
......
from enum import Enum, auto
class GpuMode(Enum):
AUTO = auto()
FULL = auto()
NONE = auto()
KERNEL = auto()
def get_best_gpu_mode(model):
""" Returns the best gpu mode with respect to the model that is to be run, gpu resources, ...
TODO
:return: GpuMode
"""
from core.models.abstract_statistical_model import AbstractStatisticalModel
from core.models.longitudinal_atlas import LongitudinalAtlas
assert isinstance(model, AbstractStatisticalModel)
gpu_mode = GpuMode.FULL
if model.name == "LongitudinalAtlas":
gpu_mode = GpuMode.KERNEL
return gpu_mode
import os
import torch
from core import GpuMode
from support import utilities
logger_format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
......@@ -63,7 +63,8 @@ scale_initial_step_size = True
downsampling_factor = 1
dense_mode = False
use_cuda = True if torch.cuda.is_available() else False
gpu_mode = GpuMode.KERNEL
# use_cuda = True if torch.cuda.is_available() else False
_cuda_is_used = False # true if at least one operation will use CUDA.
_keops_is_used = False # true if at least one keops kernel operation will take place.
......@@ -138,7 +139,7 @@ def update_dtype(new_dtype):
tensor_integer_type = utilities.get_torch_integer_type(dtype)
def update_use_cuda(new_use_cuda):
global use_cuda
assert isinstance(new_use_cuda, bool)
use_cuda = new_use_cuda
# def update_use_cuda(new_use_cuda):
# global use_cuda
# assert isinstance(new_use_cuda, bool)
# use_cuda = new_use_cuda
import numpy as np
import torch
import logging
......@@ -87,12 +86,12 @@ class MultiObjectAttachment:
Compute the current distance between source and target, assuming points are the new points of the source
We assume here that the target never moves.
"""
c1, n1, c2, n2 = MultiObjectAttachment.__get_source_and_target_centers_and_normals(points, source, target)
device, _ = utilities.get_best_device(kernel.gpu_mode)
c1, n1, c2, n2 = MultiObjectAttachment.__get_source_and_target_centers_and_normals(points, source, target, device=device)
def current_scalar_product(points_1, points_2, normals_1, normals_2):
assert points_1.device == points_2.device == normals_1.device == normals_2.device, 'tensors must be on the same device'
return torch.dot(normals_1.view(-1), kernel.convolve(points_1, points_2, normals_2, return_to_cpu=False).view(-1))
return torch.dot(normals_1.view(-1), kernel.convolve(points_1, points_2, normals_2).view(-1))
if target.norm is None:
target.norm = current_scalar_product(c2, c2, n2, n2)
......@@ -105,12 +104,12 @@ class MultiObjectAttachment:
Compute the point cloud distance between source and target, assuming points are the new points of the source
We assume here that the target never moves.
"""
c1, n1, c2, n2 = MultiObjectAttachment.__get_source_and_target_centers_and_normals(points, source, target)
device, _ = utilities.get_best_device(kernel.gpu_mode)
c1, n1, c2, n2 = MultiObjectAttachment.__get_source_and_target_centers_and_normals(points, source, target, device=device)
def point_cloud_scalar_product(points_1, points_2, normals_1, normals_2):
return torch.dot(normals_1.view(-1),
kernel.convolve(points_1, points_2, normals_2, mode='pointcloud', return_to_cpu=False).view(-1))
kernel.convolve(points_1, points_2, normals_2, mode='pointcloud').view(-1))
if target.norm is None:
target.norm = point_cloud_scalar_product(c2, c2, n2, n2)
......@@ -125,8 +124,8 @@ class MultiObjectAttachment:
source and target are SurfaceMesh objects
points are source points (torch tensor)
"""
c1, n1, c2, n2 = MultiObjectAttachment.__get_source_and_target_centers_and_normals(points, source, target)
device, _ = utilities.get_best_device(kernel.gpu_mode)
c1, n1, c2, n2 = MultiObjectAttachment.__get_source_and_target_centers_and_normals(points, source, target, device=device)
# alpha = normales non unitaires
areaa = torch.norm(n1, 2, 1)
......@@ -136,7 +135,7 @@ class MultiObjectAttachment:
nbeta = n2 / areab.unsqueeze(1)
def varifold_scalar_product(x, y, areaa, areab, nalpha, nbeta):
return torch.dot(areaa.view(-1), kernel.convolve((x, nalpha), (y, nbeta), areab.view(-1, 1), mode='varifold', return_to_cpu=False).view(-1))
return torch.dot(areaa.view(-1), kernel.convolve((x, nalpha), (y, nbeta), areab.view(-1, 1), mode='varifold').view(-1))
if target.norm is None:
target.norm = varifold_scalar_product(c2, c2, areab, areab, nbeta, nbeta)
......@@ -175,17 +174,19 @@ class MultiObjectAttachment:
####################################################################################################################
@staticmethod
def __get_source_and_target_centers_and_normals(points, source, target):
def __get_source_and_target_centers_and_normals(points, source, target, device=None):
if device is None:
device = points.device
dtype = str(points.dtype)
use_cuda = points.device.type == 'cuda'
c1, n1 = source.get_centers_and_normals(points,
tensor_scalar_type=utilities.get_torch_scalar_type(dtype=dtype),
tensor_integer_type=utilities.get_torch_integer_type(dtype=dtype),
device=points.device)
device=device)
c2, n2 = target.get_centers_and_normals(tensor_scalar_type=utilities.get_torch_scalar_type(dtype=dtype),
tensor_integer_type=utilities.get_torch_integer_type(dtype=dtype),
device=points.device)
device=device)
assert c1.device == n1.device == c2.device == n2.device, 'all tensors must be on the same device, c1.device=' + str(c1.device) \
+ ', n1.device=' + str(n1.device)\
......
import time
import warnings
from copy import deepcopy
import support.kernels as kernel_factory
......@@ -40,7 +39,7 @@ class Exponential:
self.kernel = kernel
if shoot_kernel_type is not None:
self.shoot_kernel = kernel_factory.factory(shoot_kernel_type, kernel_width=kernel.kernel_width)
self.shoot_kernel = kernel_factory.factory(shoot_kernel_type, gpu_mode=kernel.gpu_mode, kernel_width=kernel.kernel_width)
else:
self.shoot_kernel = self.kernel
......@@ -73,10 +72,14 @@ class Exponential:
# self.cholesky_matrices = {}
def move_data_to_(self, device):
self.initial_control_points = utilities.move_data(self.initial_control_points, device)
self.initial_momenta = utilities.move_data(self.initial_momenta, device)
self.initial_template_points = {key: utilities.move_data(value, device) for key, value in
self.initial_template_points.items()}
if self.initial_control_points is not None:
self.initial_control_points = utilities.move_data(self.initial_control_points, device)
if self.initial_momenta is not None:
self.initial_momenta = utilities.move_data(self.initial_momenta, device)
if self.initial_template_points is not None:
self.initial_template_points = {key: utilities.move_data(value, device) for key, value in
self.initial_template_points.items()}
def light_copy(self):
light_copy = Exponential(self.dense_mode,
......@@ -139,7 +142,7 @@ class Exponential:
"""
returns the scalar product 'mom1 K(cp) mom 2'
"""
return torch.sum(mom1 * self.kernel.convolve(cp, cp, mom2, return_to_cpu=False))
return torch.sum(mom1 * self.kernel.convolve(cp, cp, mom2))
def get_template_points(self, time_index=None):
"""
......@@ -239,8 +242,7 @@ class Exponential:
landmark_points = [self.initial_template_points['landmark_points']]
for i in range(self.number_of_time_points - 1):
d_pos = self.kernel.convolve(landmark_points[i], self.control_points_t[i], self.momenta_t[i],
return_to_cpu=False)
d_pos = self.kernel.convolve(landmark_points[i], self.control_points_t[i], self.momenta_t[i])
landmark_points.append(landmark_points[i] + dt * d_pos)
if self.use_rk2_for_flow:
......@@ -249,14 +251,12 @@ class Exponential:
if i < self.number_of_time_points - 2:
landmark_points[-1] = landmark_points[i] + dt / 2 * \
(self.kernel.convolve(landmark_points[i + 1],
self.control_points_t[i + 1], self.momenta_t[i + 1],
return_to_cpu=False) + d_pos)
self.control_points_t[i + 1], self.momenta_t[i + 1]) + d_pos)
else:
final_cp, final_mom = self._rk2_step(self.kernel, self.control_points_t[-1], self.momenta_t[-1],
dt, return_mom=True)
landmark_points[-1] = landmark_points[i] + dt / 2 * (
self.kernel.convolve(landmark_points[i + 1], final_cp, final_mom,
return_to_cpu=False) + d_pos)
self.kernel.convolve(landmark_points[i + 1], final_cp, final_mom) + d_pos)
self.template_points_t['landmark_points'] = landmark_points
......@@ -269,8 +269,7 @@ class Exponential:
for i in range(self.number_of_time_points - 1):
vf = self.kernel.convolve(image_points[0].contiguous().view(-1, dimension), self.control_points_t[i],
self.momenta_t[i],
return_to_cpu=False).view(image_shape)
self.momenta_t[i]).view(image_shape)
dY = self._compute_image_explicit_euler_step_at_order_1(image_points[i], vf)
image_points.append(image_points[i] - dt * dY)
......@@ -483,8 +482,8 @@ class Exponential:
"""
assert cp.device == mom.device, 'tensors must be on the same device, cp.device=' + str(
cp.device) + ', mom.device=' + str(mom.device)
return cp + h * kernel.convolve(cp, cp, mom, return_to_cpu=False), \
mom - h * kernel.convolve_gradient(mom, cp, return_to_cpu=False)
return cp + h * kernel.convolve(cp, cp, mom), \
mom - h * kernel.convolve_gradient(mom, cp)
@staticmethod
def _rk2_step(kernel, cp, mom, h, return_mom=True):
......@@ -496,13 +495,13 @@ class Exponential:
assert cp.device == mom.device, 'tensors must be on the same device, cp.device=' + str(
cp.device) + ', mom.device=' + str(mom.device)
mid_cp = cp + h / 2. * kernel.convolve(cp, cp, mom, return_to_cpu=False)
mid_mom = mom - h / 2. * kernel.convolve_gradient(mom, cp, return_to_cpu=False)
mid_cp = cp + h / 2. * kernel.convolve(cp, cp, mom)
mid_mom = mom - h / 2. * kernel.convolve_gradient(mom, cp)
if return_mom:
return cp + h * kernel.convolve(mid_cp, mid_cp, mid_mom, return_to_cpu=False), \
mom - h * kernel.convolve_gradient(mid_mom, mid_cp, return_to_cpu=False)
return cp + h * kernel.convolve(mid_cp, mid_cp, mid_mom), \
mom - h * kernel.convolve_gradient(mid_mom, mid_cp)
else:
return cp + h * kernel.convolve(mid_cp, mid_cp, mid_mom, return_to_cpu=False)
return cp + h * kernel.convolve(mid_cp, mid_cp, mid_mom)
# TODO. Wrap pytorch of an efficient C code ? Use keops ? Called ApplyH in PyCa. Check Numba as well.
# @jit(parallel=True)
......
import logging
import resource
import time
import warnings
import torch
from core import default
from core.model_tools.deformations.exponential import Exponential
from in_out.array_readers_and_writers import *
import torch.multiprocessing as mp
from support import utilities
import logging
logger = logging.getLogger(__name__)
......@@ -165,7 +164,7 @@ class Geodesic:
msg = "Asking for deformed template data but the geodesic was modified and not updated"
warnings.warn(msg)
times = self._get_times()
times = self.get_times()
# Deal with the special case of a geodesic reduced to a single point.
if len(times) == 1:
......@@ -190,9 +189,12 @@ class Geodesic:
# assert times[j-1] <= time
# assert times[j] >= time
weight_left = torch.Tensor([(times[j] - time) / (times[j] - times[j - 1])]).type(self.momenta_t0.type())
weight_right = torch.Tensor([(time - times[j - 1]) / (times[j] - times[j - 1])]).type(self.momenta_t0.type())
template_t = self._get_template_points_trajectory()
device, _ = utilities.get_best_device(self.backward_exponential.kernel.gpu_mode)
weight_left = utilities.move_data([(times[j] - time) / (times[j] - times[j - 1])], device=device, dtype=self.momenta_t0.dtype)
weight_right = utilities.move_data([(time - times[j - 1]) / (times[j] - times[j - 1])], device=device, dtype=self.momenta_t0.dtype)
template_t = {key: [utilities.move_data(v, device=device) for v in value] for key, value in self.get_template_points_trajectory().items()}
deformed_points = {key: weight_left * value[j - 1] + weight_right * value[j]
for key, value in template_t.items()}
......@@ -213,6 +215,8 @@ class Geodesic:
if self.shoot_is_modified or self.flow_is_modified:
device, _ = utilities.get_best_device(self.backward_exponential.kernel.gpu_mode)
# Backward exponential -------------------------------------------------------------------------------------
length = self.t0 - self.tmin
self.backward_exponential.number_of_time_points = \
......@@ -223,6 +227,7 @@ class Geodesic:
if self.flow_is_modified:
self.backward_exponential.set_initial_template_points(self.template_points_t0)
if self.backward_exponential.number_of_time_points > 1:
self.backward_exponential.move_data_to_(device=device)
self.backward_exponential.update()
# Forward exponential --------------------------------------------------------------------------------------
......@@ -235,6 +240,7 @@ class Geodesic:
if self.flow_is_modified:
self.forward_exponential.set_initial_template_points(self.template_points_t0)
if self.forward_exponential.number_of_time_points > 1:
self.forward_exponential.move_data_to_(device=device)
self.forward_exponential.update()
self.shoot_is_modified = False
......@@ -312,11 +318,7 @@ class Geodesic:
+ parallel_transport_t + parallel_transport_t_forward_extension[1:]
return parallel_transport_t
####################################################################################################################