From 3387656bdd7e8db1275882a424cc4861a5726751 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 13:22:29 +0100 Subject: [PATCH 01/70] ENH: Add (initial) support for tensorflow interconnectivity --- .../operator/tensorflow_cnn_tomography.py | 30 +++++++++++ examples/operator/tensorflow_layer_matrix.py | 32 ++++++++++++ .../operator/tensorflow_layer_tomography.py | 32 ++++++++++++ odl/operator/oputils.py | 51 ++++++++++++++++++- 4 files changed, 144 insertions(+), 1 deletion(-) create mode 100644 examples/operator/tensorflow_cnn_tomography.py create mode 100644 examples/operator/tensorflow_layer_matrix.py create mode 100644 examples/operator/tensorflow_layer_tomography.py diff --git a/examples/operator/tensorflow_cnn_tomography.py b/examples/operator/tensorflow_cnn_tomography.py new file mode 100644 index 00000000000..aecd8152f3b --- /dev/null +++ b/examples/operator/tensorflow_cnn_tomography.py @@ -0,0 +1,30 @@ +"""Example of how to convert an ODL operator to a tensorflow layer.""" + +import tensorflow as tf +import numpy as np +import odl + +with tf.Session() as sess: + tf.global_variables_initializer().run() + + space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space) + ray_transform = odl.tomo.RayTransform(space, geometry) + + x = tf.constant(np.asarray(ray_transform.domain.one())) + z = tf.constant(np.asarray(ray_transform.range.one())) + + odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') + y = odl_op_layer(x) + + # Evaluate using tensorflow + print(y.eval()) + + # Compare result with pure ODL + print(ray_transform(x.eval())) + + # Evaluate the adjoint of the derivative, called gradient in tensorflow + print(tf.gradients(y, [x], z)[0].eval()) + + # Compare result with pure ODL + print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/examples/operator/tensorflow_layer_matrix.py b/examples/operator/tensorflow_layer_matrix.py new file mode 100644 index 00000000000..7780230a072 --- /dev/null +++ b/examples/operator/tensorflow_layer_matrix.py @@ -0,0 +1,32 @@ +"""Example of how to convert an ODL operator to a tensorflow layer.""" + +import tensorflow as tf +import numpy as np +import odl + + +with tf.Session() as sess: + tf.global_variables_initializer().run() + + matrix = np.array([[1, 2], + [0, 0], + [0, 1]], dtype='float32') + odl_op = odl.MatrixOperator(matrix) + + x = tf.constant([1., 2.]) + z = tf.constant([1., 2., 3.]) + + odl_op_layer = odl.as_tensorflow_layer(odl_op, 'MatrixOperator') + y = odl_op_layer(x) + + # Evaluate using tensorflow + print(y.eval()) + + # Compare result with pure ODL + print(odl_op(x.eval())) + + # Evaluate the adjoint of the derivative, called gradient in tensorflow + print(tf.gradients(y, [x], z)[0].eval()) + + # Compare result with pure ODL + print(odl_op.derivative(x.eval()).adjoint(z.eval())) diff --git a/examples/operator/tensorflow_layer_tomography.py b/examples/operator/tensorflow_layer_tomography.py new file mode 100644 index 00000000000..66ad0e2124c --- /dev/null +++ b/examples/operator/tensorflow_layer_tomography.py @@ -0,0 +1,32 @@ +"""Example of how to convert an ODL ray transform to a tensorflow layer.""" + +import tensorflow as tf +import numpy as np +import odl + + +with tf.Session() as sess: + tf.global_variables_initializer().run() + + space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], + dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space) + ray_transform = odl.tomo.RayTransform(space, geometry) + + x = tf.constant(np.asarray(ray_transform.domain.one())) + z = tf.constant(np.asarray(ray_transform.range.one())) + + odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') + y = odl_op_layer(x) + + # Evaluate using tensorflow + print(y.eval()) + + # Compare result with pure ODL + print(ray_transform(x.eval())) + + # Evaluate the adjoint of the derivative, called gradient in tensorflow + print(tf.gradients(y, [x], z)[0].eval()) + + # Compare result with pure ODL + print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index aeabd792048..a8ace4559b8 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -21,7 +21,8 @@ from odl.util import as_flat_array __all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator', - 'as_scipy_functional', 'as_proximal_lang_operator') + 'as_scipy_functional', 'as_proximal_lang_operator', + 'as_tensorflow_layer') def matrix_representation(op): @@ -439,6 +440,54 @@ def adjoint(inp, out): adjoint=adjoint, norm_bound=norm_bound) + +def as_tensorflow_layer(odl_op, name): + """Convert ``Operator`` to tensorflow layer.""" + import tensorflow as tf + from tensorflow.python.framework import ops + + def py_func(func, inp, Tout, stateful=True, name=None, grad=None): + """Define custom py_func which takes also a grad op as argument.""" + if grad is None: + return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + else: + # Need to generate a unique name to avoid duplicates: + rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8)) + + tf.RegisterGradient(rnd_name)(grad) + g = tf.get_default_graph() + with g.gradient_override_map({"PyFunc": rnd_name}): + return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + + def tensorflow_layer_grad_impl(x, dx): + def _impl(x, dx): + return np.asarray(odl_op.derivative(x).adjoint(dx)) + + with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): + result = py_func(_impl, + [x, dx], + [tf.float32]) + return result[0] + + # Actual gradient: + def tensorflow_layer_grad(op, grad): + x = op.inputs[0] + return tensorflow_layer_grad_impl(x, grad) + + def tensorflow_layer(x): + def _impl(x): + return np.asarray(odl_op(x)) + + with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: + result = py_func(_impl, + [x], + [tf.float32], + name=name_call, + grad=tensorflow_layer_grad) + return result[0] + + return tensorflow_layer + if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From 3ab822b038a9cde49e351014332ad0b276a652a0 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 13:40:30 +0100 Subject: [PATCH 02/70] ENH: Add more tensorflow examples --- .../operator/tensorflow_cnn_tomography.py | 30 -------------- ...y.py => tensorflow_layer_ray_transform.py} | 0 examples/operator/tensorflow_tomography.py | 39 +++++++++++++++++++ 3 files changed, 39 insertions(+), 30 deletions(-) delete mode 100644 examples/operator/tensorflow_cnn_tomography.py rename examples/operator/{tensorflow_layer_tomography.py => tensorflow_layer_ray_transform.py} (100%) create mode 100644 examples/operator/tensorflow_tomography.py diff --git a/examples/operator/tensorflow_cnn_tomography.py b/examples/operator/tensorflow_cnn_tomography.py deleted file mode 100644 index aecd8152f3b..00000000000 --- a/examples/operator/tensorflow_cnn_tomography.py +++ /dev/null @@ -1,30 +0,0 @@ -"""Example of how to convert an ODL operator to a tensorflow layer.""" - -import tensorflow as tf -import numpy as np -import odl - -with tf.Session() as sess: - tf.global_variables_initializer().run() - - space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space) - ray_transform = odl.tomo.RayTransform(space, geometry) - - x = tf.constant(np.asarray(ray_transform.domain.one())) - z = tf.constant(np.asarray(ray_transform.range.one())) - - odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') - y = odl_op_layer(x) - - # Evaluate using tensorflow - print(y.eval()) - - # Compare result with pure ODL - print(ray_transform(x.eval())) - - # Evaluate the adjoint of the derivative, called gradient in tensorflow - print(tf.gradients(y, [x], z)[0].eval()) - - # Compare result with pure ODL - print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/examples/operator/tensorflow_layer_tomography.py b/examples/operator/tensorflow_layer_ray_transform.py similarity index 100% rename from examples/operator/tensorflow_layer_tomography.py rename to examples/operator/tensorflow_layer_ray_transform.py diff --git a/examples/operator/tensorflow_tomography.py b/examples/operator/tensorflow_tomography.py new file mode 100644 index 00000000000..bea55677564 --- /dev/null +++ b/examples/operator/tensorflow_tomography.py @@ -0,0 +1,39 @@ +"""Example of how to solve a simple tomography problem using tensorflow.""" + +import tensorflow as tf +import numpy as np +import odl + + +with tf.Session() as sess: + # Create ODL data structures + space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], + dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space) + ray_transform = odl.tomo.RayTransform(space, geometry) + phantom = odl.phantom.shepp_logan(space, True) + data = ray_transform(phantom) + + # Create tensorflow layer from odl operator + odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') + x = tf.Variable(tf.constant(0.0, shape=space.shape), name="x") + + # Create constant right hand side + y = tf.constant(np.asarray(data)) + + # Define l2 loss function + loss = tf.nn.l2_loss(odl_op_layer(x) - y) + + # Train using the adam optimizer + optimizer = tf.train.AdamOptimizer(1e-1).minimize(loss) + + # Initialize all TF variables + tf.global_variables_initializer().run() + + # Solve with an ODL callback to see what happens + callback = odl.solvers.CallbackShow() + + for i in range(100): + sess.run(optimizer) + + callback(space.element(x.eval())) From 83c27e40797bc0b5ad6fb00132e4d5b8f39fc901 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 17:27:01 +0100 Subject: [PATCH 03/70] ENH: Improve handling of shapes of tensorflow layers --- odl/operator/oputils.py | 82 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 74 insertions(+), 8 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index a8ace4559b8..86d22371e64 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -441,8 +441,29 @@ def adjoint(inp, out): norm_bound=norm_bound) -def as_tensorflow_layer(odl_op, name): - """Convert ``Operator`` to tensorflow layer.""" +def as_tensorflow_layer(odl_op, name='ODLOperator'): + """Convert ``Operator`` to tensorflow layer. + + Parameters + ---------- + odl_op : `Operator` + The operator that should be wrapped to a tensorflow layer. + name : str + Tensorflow name of the operator + + Returns + ------- + tensorflow_layer : callable + Callable that, when called with an `tensorflow.Tensor` of shape + `(n, *odl_op.domain.shape, 1)` returns a tensor of shape + `(n, *odl_op.range.shape, 1)` where ``n`` is the number of batches. + Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. + The `dtype` of the tensor is the same as the respective ODL spaces. + + If ``odl_op`` implements `Operator.derivative` which in turn implements + `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and + gradients propagate as expected. + """ import tensorflow as tf from tensorflow.python.framework import ops @@ -457,26 +478,66 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): tf.RegisterGradient(rnd_name)(grad) g = tf.get_default_graph() with g.gradient_override_map({"PyFunc": rnd_name}): - return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + return tf.py_func(func, inp, Tout, stateful=stateful, + name=name) def tensorflow_layer_grad_impl(x, dx): + """Implementation of the tensorflow gradient. + + Gradient in tensorflow is equivalent to the adjoint of the derivative + in ODL. + """ + x_shape = x.get_shape() + dx_shape = dx.get_shape() + out_shape = (int(x_shape[0]),) + odl_op.domain.shape + (1,) + + # Validate input shape + assert x_shape[0] == dx_shape[0] + assert x_shape[1:] == odl_op.domain.shape + (1,) + assert dx_shape[1:] == odl_op.range.shape + (1,) + def _impl(x, dx): - return np.asarray(odl_op.derivative(x).adjoint(dx)) + out = np.empty(out_shape, odl_op.domain.dtype) + for i in range(x_shape[0]): + xi = x[i, ..., 0] + dxi = dx[i, ..., 0] + result = odl_op.derivative(xi).adjoint(dxi) + out[i] = np.asarray(result)[..., None] + return out with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): result = py_func(_impl, [x, dx], [tf.float32]) - return result[0] - # Actual gradient: + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result + def tensorflow_layer_grad(op, grad): + """Thin wrapper for the gradient.""" x = op.inputs[0] return tensorflow_layer_grad_impl(x, grad) def tensorflow_layer(x): + """Implementation of the tensorflow call. + + Gradient in tensorflow is equivalent to the adjoint of the derivative + in ODL. + """ + x_shape = x.get_shape() + out_shape = (int(x_shape[0]),) + odl_op.range.shape + (1,) + + # Validate input shape + assert x_shape[1:] == odl_op.domain.shape + (1,) + def _impl(x): - return np.asarray(odl_op(x)) + out = np.empty(out_shape, odl_op.range.dtype) + for i in range(x_shape[0]): + out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] + return out with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: result = py_func(_impl, @@ -484,7 +545,12 @@ def _impl(x): [tf.float32], name=name_call, grad=tensorflow_layer_grad) - return result[0] + + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result return tensorflow_layer From b4b4a2077e8b53b6c1bf4c8eb5d7cfb16934f02c Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 17:27:11 +0100 Subject: [PATCH 04/70] ENH: Add more tensorflow examples --- .../operator/tensorflow_tomography_cnn.py | 61 +++++++++++++ .../tensorflow_tomography_cnn_gradient.py | 91 +++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 examples/operator/tensorflow_tomography_cnn.py create mode 100644 examples/operator/tensorflow_tomography_cnn_gradient.py diff --git a/examples/operator/tensorflow_tomography_cnn.py b/examples/operator/tensorflow_tomography_cnn.py new file mode 100644 index 00000000000..2e52ef29d9d --- /dev/null +++ b/examples/operator/tensorflow_tomography_cnn.py @@ -0,0 +1,61 @@ +"""Example of how to solve a simple tomography problem using tensorflow.""" + +import tensorflow as tf +import numpy as np +import odl + + +def conv2d(x, W): + return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME') + + +with tf.Session() as sess: + # Create ODL data structures + space = odl.uniform_discr([-64, -64], [64, 64], [512, 512], + dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space, angles=20) + ray_transform = odl.tomo.RayTransform(space, geometry) + phantom = odl.phantom.shepp_logan(space, True) + data = ray_transform(phantom) + fbp = odl.tomo.fbp_op(ray_transform)(data) + + # Create tensorflow layer from odl operator + odl_op_layer = odl.as_tensorflow_layer(ray_transform, + 'RayTransform') + x = tf.constant(np.asarray(fbp), name="x") + x = x[None, :, :, None] + + # Create constant right hand side + y = tf.constant(np.asarray(data)) + + init = np.random.randn(*[3, 3, 1, 32]) * 0.1 + W = tf.Variable(tf.constant(init, dtype='float32')) + x = tf.nn.relu(conv2d(x, W)) + + for i in range(2): + init = np.random.randn(*[3, 3, 32, 32]) * 0.1 + W = tf.Variable(tf.constant(init, dtype='float32')) + x = tf.nn.relu(conv2d(x, W)) + + init = np.random.randn(*[3, 3, 32, 1]) * 0.1 + W = tf.Variable(tf.constant(init, dtype='float32')) + x = tf.nn.relu(conv2d(x, W)) + + # Reshape for ODL + x = x[0, :, :, 0] + + # Define l2 loss function + loss = tf.nn.l2_loss(odl_op_layer(x) - y) + + # Train using the adam optimizer + optimizer = tf.train.AdamOptimizer(1e-3).minimize(loss) + + # Initialize all TF variables + tf.global_variables_initializer().run() + + # Solve with an ODL callback to see what happens + callback = odl.solvers.CallbackShow(clim=[0, 1]) + + for i in range(1000): + sess.run(optimizer) + callback(space.element(x.eval())) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py new file mode 100644 index 00000000000..e2636a64caf --- /dev/null +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -0,0 +1,91 @@ +"""Example of how to solve a simple tomography problem using tensorflow.""" + +import tensorflow as tf +import numpy as np +import odl + + +def conv2d(x, W): + return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME') + + +def var(x): + return tf.Variable(tf.constant(x, dtype='float32')) + + +with tf.Session() as sess: + # Create ODL data structures + space = odl.uniform_discr([-64, -64], [64, 64], [256, 256], + dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space) + ray_transform = odl.tomo.RayTransform(space, geometry) + phantom = odl.phantom.shepp_logan(space, True) + data = ray_transform(phantom) + noisy_data = data + odl.phantom.white_noise(ray_transform.range) * np.mean(data) * 0.1 + fbp = odl.tomo.fbp_op(ray_transform)(noisy_data) + + # Create tensorflow layer from odl operator + odl_op_layer = odl.as_tensorflow_layer(ray_transform, + 'RayTransform') + odl_op_layer_adjoint = odl.as_tensorflow_layer(ray_transform.adjoint, + 'RayTransformAdjoint') + x = tf.constant(np.asarray(fbp)[None, :, :, None], name="x") + + x_true = tf.constant(np.asarray(phantom)[None, :, :, None], name="x_true") + + s = tf.Variable(tf.constant(np.zeros([1, 256, 256, 5]), dtype='float32')) + + # Create constant right hand side + y = tf.constant(np.asarray(noisy_data)[None, :, :, None]) + + w1 = var(np.random.randn(*[3, 3, 7, 32]) * 0.01) + b1 = var(np.random.randn(*[1, 256, 256, 32]) * 0.001) + + w2 = var(np.random.randn(*[3, 3, 32, 32]) * 0.01) + b2 = var(np.random.randn(*[1, 256, 256, 32]) * 0.001) + + w3 = var(np.random.randn(*[3, 3, 32, 6]) * 0.01) + b3 = var(np.random.randn(*[1, 256, 256, 6]) * 0.001) + + n_iter = 5 + x_vals = [] + + loss = 0 + for i in range(n_iter): + # Reshape for ODL + residual = odl_op_layer(x) - y + dx = odl_op_layer_adjoint(residual) + + dx = tf.concat([x, dx, s], axis=3) + + dx = tf.nn.relu(conv2d(dx, w1) + b1) + dx = tf.nn.relu(conv2d(dx, w2) + b2) + dx = tf.nn.dropout(dx, 0.95) + dx = tf.nn.tanh(conv2d(dx, w3) + b3) + + s = dx[..., 1:] + dx = dx[..., 0][..., None] + + x = x + dx + + x_vals.append(x) + + loss = loss + tf.nn.l2_loss(x - x_true) * (2. ** (i - n_iter)) + + # Train using the adam optimizer + learning_rate = tf.placeholder(tf.float32, shape=[]) + optimizer = tf.train.AdamOptimizer(learning_rate).minimize(loss) + + # Initialize all TF variables + tf.global_variables_initializer().run() + + # Solve with an ODL callback to see what happens + callback = odl.solvers.CallbackShow(clim=[0, 1]) + + for i in range(1000): + if i < 100: + sess.run(optimizer, feed_dict={learning_rate: 0.001}) + else: + sess.run(optimizer, feed_dict={learning_rate: 0.0001}) + callback((space ** n_iter).element([xi.eval() for xi in x_vals])) + print(loss.eval()) From 8d500aa2c8282e4f863436a1de7f5233c3082df2 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 21:39:50 +0100 Subject: [PATCH 05/70] ENH: Optimizations to the tensorflow interface --- odl/operator/oputils.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 86d22371e64..0caf6b2a026 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -503,12 +503,17 @@ def _impl(x, dx): dxi = dx[i, ..., 0] result = odl_op.derivative(xi).adjoint(dxi) out[i] = np.asarray(result)[..., None] + + # TODO: Improve + scale = odl_op.domain.cell_volume / odl_op.range.cell_volume + out *= scale return out with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): result = py_func(_impl, [x, dx], - [tf.float32]) + [tf.float32], + stateful=False) # We must manually set the output shape since tensorflow cannot # figure it out @@ -544,6 +549,7 @@ def _impl(x): [x], [tf.float32], name=name_call, + stateful=False, grad=tensorflow_layer_grad) # We must manually set the output shape since tensorflow cannot From faef60fbc0d937d583e95efac5cac3451a44fab0 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 21:40:14 +0100 Subject: [PATCH 06/70] ENH: Improved CNN gradient tomography --- .../tensorflow_tomography_cnn_gradient.py | 89 ++++++++++++------- 1 file changed, 58 insertions(+), 31 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index e2636a64caf..bd9bca97733 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -5,6 +5,19 @@ import odl +def random_ellipse(): + return (np.random.rand() - 0.3, + np.random.exponential() * 0.2, np.random.exponential() * 0.2, + np.random.rand() - 0.5, np.random.rand() - 0.5, + np.random.rand() * 2 * np.pi) + + +def random_phantom(spc): + n = np.random.poisson(10) + ellipses = [random_ellipse() for _ in range(n)] + return odl.phantom.ellipsoid_phantom(spc, ellipses) + + def conv2d(x, W): return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME') @@ -15,62 +28,74 @@ def var(x): with tf.Session() as sess: # Create ODL data structures - space = odl.uniform_discr([-64, -64], [64, 64], [256, 256], + size = 128 + space = odl.uniform_discr([-64, -64], [64, 64], [size, size], dtype='float32') geometry = odl.tomo.parallel_beam_geometry(space) ray_transform = odl.tomo.RayTransform(space, geometry) - phantom = odl.phantom.shepp_logan(space, True) - data = ray_transform(phantom) - noisy_data = data + odl.phantom.white_noise(ray_transform.range) * np.mean(data) * 0.1 - fbp = odl.tomo.fbp_op(ray_transform)(noisy_data) + # Create tensorflow layer from odl operator odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') odl_op_layer_adjoint = odl.as_tensorflow_layer(ray_transform.adjoint, 'RayTransformAdjoint') - x = tf.constant(np.asarray(fbp)[None, :, :, None], name="x") - x_true = tf.constant(np.asarray(phantom)[None, :, :, None], name="x_true") + n_data = 50 + x_arr = np.empty((n_data, space.shape[0], space.shape[1], 1), dtype='float32') + y_arr = np.empty((n_data, ray_transform.range.shape[0], ray_transform.range.shape[1], 1), dtype='float32') + x_true_arr = np.empty((n_data, space.shape[0], space.shape[1], 1), dtype='float32') + + for i in range(n_data): + if i == 0: + phantom = odl.phantom.shepp_logan(space, True) + else: + phantom = random_phantom(space) + data = ray_transform(phantom) + noisy_data = data + odl.phantom.white_noise(ray_transform.range) * np.mean(data) * 0.1 + fbp = odl.tomo.fbp_op(ray_transform)(noisy_data) + + x_arr[i] = np.asarray(fbp)[..., None] + x_true_arr[i] = np.asarray(phantom)[..., None] + y_arr[i] = np.asarray(noisy_data)[..., None] + + x = tf.constant(x_arr, name="x") + x_true = tf.constant(x_true_arr, name="x_true") + y = tf.constant(y_arr, name='y') - s = tf.Variable(tf.constant(np.zeros([1, 256, 256, 5]), dtype='float32')) + s = tf.Variable(tf.constant(np.zeros([n_data, size, size, 5]), dtype='float32')) # Create constant right hand side - y = tf.constant(np.asarray(noisy_data)[None, :, :, None]) - w1 = var(np.random.randn(*[3, 3, 7, 32]) * 0.01) - b1 = var(np.random.randn(*[1, 256, 256, 32]) * 0.001) + w1 = tf.Variable(tf.truncated_normal([3, 3, 7, 32], stddev=0.01)) + b1 = var(np.ones(32) * 0.1) - w2 = var(np.random.randn(*[3, 3, 32, 32]) * 0.01) - b2 = var(np.random.randn(*[1, 256, 256, 32]) * 0.001) + w2 = tf.Variable(tf.truncated_normal([3, 3, 32, 32], stddev=0.01)) + b2 = var(np.ones(32) * 0.1) - w3 = var(np.random.randn(*[3, 3, 32, 6]) * 0.01) - b3 = var(np.random.randn(*[1, 256, 256, 6]) * 0.001) + w3 = tf.Variable(tf.truncated_normal([3, 3, 32, 6], stddev=0.01)) + b3 = var(np.random.randn(6) * 0.001) n_iter = 5 - x_vals = [] loss = 0 for i in range(n_iter): # Reshape for ODL - residual = odl_op_layer(x) - y - dx = odl_op_layer_adjoint(residual) + gradx = odl_op_layer_adjoint(odl_op_layer(x) - y) - dx = tf.concat([x, dx, s], axis=3) + dx = tf.concat([x, gradx, s], axis=3) dx = tf.nn.relu(conv2d(dx, w1) + b1) dx = tf.nn.relu(conv2d(dx, w2) + b2) - dx = tf.nn.dropout(dx, 0.95) - dx = tf.nn.tanh(conv2d(dx, w3) + b3) + dx = tf.nn.dropout(dx, 0.8) + dx = conv2d(dx, w3) + b3 - s = dx[..., 1:] - dx = dx[..., 0][..., None] + s = tf.nn.relu(dx[..., 1:]) + dx = tf.nn.tanh(dx[..., 0][..., None]) x = x + dx - x_vals.append(x) - - loss = loss + tf.nn.l2_loss(x - x_true) * (2. ** (i - n_iter)) + loss = loss + tf.nn.l2_loss(x - x_true) * (2. ** (i - n_iter)) / n_data # Train using the adam optimizer learning_rate = tf.placeholder(tf.float32, shape=[]) @@ -80,12 +105,14 @@ def var(x): tf.global_variables_initializer().run() # Solve with an ODL callback to see what happens - callback = odl.solvers.CallbackShow(clim=[0, 1]) + callback = odl.solvers.CallbackShow() - for i in range(1000): + for i in range(10000): if i < 100: sess.run(optimizer, feed_dict={learning_rate: 0.001}) else: - sess.run(optimizer, feed_dict={learning_rate: 0.0001}) - callback((space ** n_iter).element([xi.eval() for xi in x_vals])) - print(loss.eval()) + sess.run(optimizer, feed_dict={learning_rate: 0.001}) + + if i % 10 == 0: + callback(space.element(x[0].eval())) + print(loss.eval()) From c28586fef98fca5ff824c86ed2844e381134f297 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Feb 2017 13:49:56 +0100 Subject: [PATCH 07/70] ENH: Support variable size input for tensorflow --- odl/operator/oputils.py | 42 +++++++++++++++++++++++++++++++++++------ 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 0caf6b2a026..c7073c11c3d 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -489,7 +489,15 @@ def tensorflow_layer_grad_impl(x, dx): """ x_shape = x.get_shape() dx_shape = dx.get_shape() - out_shape = (int(x_shape[0]),) + odl_op.domain.shape + (1,) + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.range.shape + (1,) + out_shape = (n_x,) + odl_op.domain.shape + (1,) # Validate input shape assert x_shape[0] == dx_shape[0] @@ -497,8 +505,15 @@ def tensorflow_layer_grad_impl(x, dx): assert dx_shape[1:] == odl_op.range.shape + (1,) def _impl(x, dx): - out = np.empty(out_shape, odl_op.domain.dtype) - for i in range(x_shape[0]): + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.domain.dtype) + for i in range(x_out_shape[0]): xi = x[i, ..., 0] dxi = dx[i, ..., 0] result = odl_op.derivative(xi).adjoint(dxi) @@ -533,14 +548,29 @@ def tensorflow_layer(x): in ODL. """ x_shape = x.get_shape() - out_shape = (int(x_shape[0]),) + odl_op.range.shape + (1,) + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.domain.shape + (1,) + out_shape = (n_x,) + odl_op.range.shape + (1,) # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) def _impl(x): - out = np.empty(out_shape, odl_op.range.dtype) - for i in range(x_shape[0]): + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.range.dtype) + for i in range(x_out_shape[0]): out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] return out From 8b0efdd98b915ec2520989a147d2cb6d1deb852f Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Feb 2017 14:27:05 +0100 Subject: [PATCH 08/70] ENH: Improved cnn gradient example --- .../tensorflow_tomography_cnn_gradient.py | 47 ++++++++++++------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index bd9bca97733..8c8578e900d 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -26,7 +26,13 @@ def var(x): return tf.Variable(tf.constant(x, dtype='float32')) +def create_variable(name, shape): + variable = tf.Variable(tf.truncated_normal(shape, stddev=0.01), name=name) + return variable + + with tf.Session() as sess: + # Create ODL data structures size = 128 space = odl.uniform_discr([-64, -64], [64, 64], [size, size], @@ -59,31 +65,32 @@ def var(x): x_true_arr[i] = np.asarray(phantom)[..., None] y_arr[i] = np.asarray(noisy_data)[..., None] - x = tf.constant(x_arr, name="x") - x_true = tf.constant(x_true_arr, name="x_true") - y = tf.constant(y_arr, name='y') + x_0 = tf.placeholder(tf.float32, shape=[None, space.shape[0], space.shape[1], 1], name="x_0") + x_true = tf.placeholder(tf.float32, shape=[None, space.shape[0], space.shape[1], 1], name="x_true") + y = tf.placeholder(tf.float32, shape=[None, ray_transform.range.shape[0], ray_transform.range.shape[1], 1], name="y") - s = tf.Variable(tf.constant(np.zeros([n_data, size, size, 5]), dtype='float32')) + s = tf.fill([tf.shape(x_0)[0], space.shape[0], space.shape[1], 5], np.float32(0.0), name="s") # Create constant right hand side - w1 = tf.Variable(tf.truncated_normal([3, 3, 7, 32], stddev=0.01)) + w1 = create_variable('w1', shape=[3, 3, 7, 32]) b1 = var(np.ones(32) * 0.1) - w2 = tf.Variable(tf.truncated_normal([3, 3, 32, 32], stddev=0.01)) + w2 = create_variable('w2', shape=[3, 3, 32, 32]) b2 = var(np.ones(32) * 0.1) - w3 = tf.Variable(tf.truncated_normal([3, 3, 32, 6], stddev=0.01)) + w3 = create_variable('w3', shape=[3, 3, 32, 6]) b3 = var(np.random.randn(6) * 0.001) n_iter = 5 + x = x_0 loss = 0 for i in range(n_iter): # Reshape for ODL gradx = odl_op_layer_adjoint(odl_op_layer(x) - y) - dx = tf.concat([x, gradx, s], axis=3) + dx = tf.concat([x, gradx, s], axis=-1) dx = tf.nn.relu(conv2d(dx, w1) + b1) dx = tf.nn.relu(conv2d(dx, w2) + b2) @@ -107,12 +114,18 @@ def var(x): # Solve with an ODL callback to see what happens callback = odl.solvers.CallbackShow() - for i in range(10000): - if i < 100: - sess.run(optimizer, feed_dict={learning_rate: 0.001}) - else: - sess.run(optimizer, feed_dict={learning_rate: 0.001}) - - if i % 10 == 0: - callback(space.element(x[0].eval())) - print(loss.eval()) + for i in range(1000): + _, loss_training = sess.run([optimizer, loss], + feed_dict={learning_rate: 0.001, + x_0: x_arr[1:], + x_true: x_true_arr[1:], + y: y_arr[1:]}) + + x_value, loss_value = sess.run([x, loss], + feed_dict={learning_rate: 0.001, + x_0: x_arr[0:1], + x_true: x_true_arr[0:1], + y: y_arr[0:1]}) + + print('iter={}, training loss={}, validation loss={}'.format(i, loss_training, loss_value)) + callback(space.element(x_value[0])) From 428087c5995fcd1dbd4669813301ef68c40a26ce Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Feb 2017 15:38:56 +0100 Subject: [PATCH 09/70] ENH: Improvements to the CNN tensorflow example --- .../tensorflow_tomography_cnn_gradient.py | 37 +++++++++---------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index 8c8578e900d..dcc314f1e2b 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -26,8 +26,8 @@ def var(x): return tf.Variable(tf.constant(x, dtype='float32')) -def create_variable(name, shape): - variable = tf.Variable(tf.truncated_normal(shape, stddev=0.01), name=name) +def create_variable(name, shape, stddev=0.01): + variable = tf.Variable(tf.truncated_normal(shape, stddev=stddev), name=name) return variable @@ -65,21 +65,21 @@ def create_variable(name, shape): x_true_arr[i] = np.asarray(phantom)[..., None] y_arr[i] = np.asarray(noisy_data)[..., None] - x_0 = tf.placeholder(tf.float32, shape=[None, space.shape[0], space.shape[1], 1], name="x_0") - x_true = tf.placeholder(tf.float32, shape=[None, space.shape[0], space.shape[1], 1], name="x_true") + x_0 = tf.placeholder(tf.float32, shape=[None, size, size, 1], name="x_0") + x_true = tf.placeholder(tf.float32, shape=[None, size, size, 1], name="x_true") y = tf.placeholder(tf.float32, shape=[None, ray_transform.range.shape[0], ray_transform.range.shape[1], 1], name="y") - s = tf.fill([tf.shape(x_0)[0], space.shape[0], space.shape[1], 5], np.float32(0.0), name="s") + s = tf.fill([tf.shape(x_0)[0], size, size, 5], np.float32(0.0), name="s") # Create constant right hand side - w1 = create_variable('w1', shape=[3, 3, 7, 32]) + w1 = tf.Variable(tf.truncated_normal([3, 3, 7, 32], stddev=0.01)) b1 = var(np.ones(32) * 0.1) - w2 = create_variable('w2', shape=[3, 3, 32, 32]) + w2 = tf.Variable(tf.truncated_normal([3, 3, 32, 32], stddev=0.01)) b2 = var(np.ones(32) * 0.1) - w3 = create_variable('w3', shape=[3, 3, 32, 6]) + w3 = tf.Variable(tf.truncated_normal([3, 3, 32, 6], stddev=0.01)) b3 = var(np.random.randn(6) * 0.001) n_iter = 5 @@ -87,19 +87,19 @@ def create_variable(name, shape): x = x_0 loss = 0 for i in range(n_iter): - # Reshape for ODL gradx = odl_op_layer_adjoint(odl_op_layer(x) - y) - dx = tf.concat([x, gradx, s], axis=-1) + update = tf.concat([x, gradx, s], axis=3) - dx = tf.nn.relu(conv2d(dx, w1) + b1) - dx = tf.nn.relu(conv2d(dx, w2) + b2) - dx = tf.nn.dropout(dx, 0.8) - dx = conv2d(dx, w3) + b3 + update = tf.nn.relu(conv2d(update, w1) + b1) + update = tf.nn.relu(conv2d(update, w2) + b2) + update = tf.nn.dropout(update, 0.8) + update = tf.nn.tanh(conv2d(update, w3) + b3) - s = tf.nn.relu(dx[..., 1:]) - dx = tf.nn.tanh(dx[..., 0][..., None]) + ds = update[..., 1:] + dx = update[..., 0][..., None] + s = s + ds x = x + dx loss = loss + tf.nn.l2_loss(x - x_true) * (2. ** (i - n_iter)) / n_data @@ -114,7 +114,7 @@ def create_variable(name, shape): # Solve with an ODL callback to see what happens callback = odl.solvers.CallbackShow() - for i in range(1000): + for i in range(10000): _, loss_training = sess.run([optimizer, loss], feed_dict={learning_rate: 0.001, x_0: x_arr[1:], @@ -122,8 +122,7 @@ def create_variable(name, shape): y: y_arr[1:]}) x_value, loss_value = sess.run([x, loss], - feed_dict={learning_rate: 0.001, - x_0: x_arr[0:1], + feed_dict={x_0: x_arr[0:1], x_true: x_true_arr[0:1], y: y_arr[0:1]}) From 787eb9974b7138d0f89aa7d12f40dd76fd03bbb8 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 3 Mar 2017 17:12:27 +0100 Subject: [PATCH 10/70] BUG: Fix bug with derivative of tensorflow operators --- odl/operator/oputils.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index c7073c11c3d..0edc9701c3f 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -477,7 +477,13 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): tf.RegisterGradient(rnd_name)(grad) g = tf.get_default_graph() - with g.gradient_override_map({"PyFunc": rnd_name}): + + if stateful: + override_name = 'PyFunc' + else: + override_name = 'PyFuncStateless' + + with g.gradient_override_map({override_name: rnd_name}): return tf.py_func(func, inp, Tout, stateful=stateful, name=name) @@ -500,7 +506,6 @@ def tensorflow_layer_grad_impl(x, dx): out_shape = (n_x,) + odl_op.domain.shape + (1,) # Validate input shape - assert x_shape[0] == dx_shape[0] assert x_shape[1:] == odl_op.domain.shape + (1,) assert dx_shape[1:] == odl_op.range.shape + (1,) @@ -508,9 +513,11 @@ def _impl(x, dx): if fixed_size: x_out_shape = out_shape assert x.shape == in_shape + assert dx.shape == out_shape else: x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == in_shape[1:] + assert x.shape[1:] == out_shape[1:] + assert dx.shape[1:] == in_shape[1:] out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): @@ -524,7 +531,7 @@ def _impl(x, dx): out *= scale return out - with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): + with ops.name_scope(name + 'Grad', "ODLTensorflowLayerGrad", [x, dx]): result = py_func(_impl, [x, dx], [tf.float32], From 4ee1f9be06ef3956d6ef2eb1b8387b5f33240daf Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 13 Mar 2017 17:41:42 +0100 Subject: [PATCH 11/70] ENH: Allow non-differentialbe operators in tensorflow --- odl/operator/oputils.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 0edc9701c3f..e6878a8a948 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -441,7 +441,7 @@ def adjoint(inp, out): norm_bound=norm_bound) -def as_tensorflow_layer(odl_op, name='ODLOperator'): +def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): """Convert ``Operator`` to tensorflow layer. Parameters @@ -450,6 +450,9 @@ def as_tensorflow_layer(odl_op, name='ODLOperator'): The operator that should be wrapped to a tensorflow layer. name : str Tensorflow name of the operator + differentiable : boolean + True if the operator should be differentiable, otherwise assumes that + the derivative is everywhere zero. Returns ------- @@ -581,13 +584,18 @@ def _impl(x): out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] return out + if differentiable: + grad = tensorflow_layer_grad + else: + grad = None + with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: result = py_func(_impl, [x], [tf.float32], name=name_call, stateful=False, - grad=tensorflow_layer_grad) + grad=grad) # We must manually set the output shape since tensorflow cannot # figure it out From 7a3a86bbe8c3be6e20cc9ee51b287e0386034e6a Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 16 Mar 2017 15:00:57 +0100 Subject: [PATCH 12/70] ENH: Add odl.contrib for less-than-optimal code --- odl/__init__.py | 1 - odl/operator/oputils.py | 167 +--------------------------------------- 2 files changed, 1 insertion(+), 167 deletions(-) diff --git a/odl/__init__.py b/odl/__init__.py index 5e3db9e4f5f..f1249139fc0 100644 --- a/odl/__init__.py +++ b/odl/__init__.py @@ -45,6 +45,5 @@ from . import datasets from . import contrib - from .util import test __all__ += ('test',) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index e6878a8a948..a16cde45e09 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -21,8 +21,7 @@ from odl.util import as_flat_array __all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator', - 'as_scipy_functional', 'as_proximal_lang_operator', - 'as_tensorflow_layer') + 'as_scipy_functional', 'as_proximal_lang_operator') def matrix_representation(op): @@ -441,170 +440,6 @@ def adjoint(inp, out): norm_bound=norm_bound) -def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): - """Convert ``Operator`` to tensorflow layer. - - Parameters - ---------- - odl_op : `Operator` - The operator that should be wrapped to a tensorflow layer. - name : str - Tensorflow name of the operator - differentiable : boolean - True if the operator should be differentiable, otherwise assumes that - the derivative is everywhere zero. - - Returns - ------- - tensorflow_layer : callable - Callable that, when called with an `tensorflow.Tensor` of shape - `(n, *odl_op.domain.shape, 1)` returns a tensor of shape - `(n, *odl_op.range.shape, 1)` where ``n`` is the number of batches. - Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. - The `dtype` of the tensor is the same as the respective ODL spaces. - - If ``odl_op`` implements `Operator.derivative` which in turn implements - `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and - gradients propagate as expected. - """ - import tensorflow as tf - from tensorflow.python.framework import ops - - def py_func(func, inp, Tout, stateful=True, name=None, grad=None): - """Define custom py_func which takes also a grad op as argument.""" - if grad is None: - return tf.py_func(func, inp, Tout, stateful=stateful, name=name) - else: - # Need to generate a unique name to avoid duplicates: - rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8)) - - tf.RegisterGradient(rnd_name)(grad) - g = tf.get_default_graph() - - if stateful: - override_name = 'PyFunc' - else: - override_name = 'PyFuncStateless' - - with g.gradient_override_map({override_name: rnd_name}): - return tf.py_func(func, inp, Tout, stateful=stateful, - name=name) - - def tensorflow_layer_grad_impl(x, dx): - """Implementation of the tensorflow gradient. - - Gradient in tensorflow is equivalent to the adjoint of the derivative - in ODL. - """ - x_shape = x.get_shape() - dx_shape = dx.get_shape() - try: - n_x = int(x_shape[0]) - fixed_size = True - except TypeError: - n_x = x_shape[0] - fixed_size = False - - in_shape = (n_x,) + odl_op.range.shape + (1,) - out_shape = (n_x,) + odl_op.domain.shape + (1,) - - # Validate input shape - assert x_shape[1:] == odl_op.domain.shape + (1,) - assert dx_shape[1:] == odl_op.range.shape + (1,) - - def _impl(x, dx): - if fixed_size: - x_out_shape = out_shape - assert x.shape == in_shape - assert dx.shape == out_shape - else: - x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == out_shape[1:] - assert dx.shape[1:] == in_shape[1:] - - out = np.empty(x_out_shape, odl_op.domain.dtype) - for i in range(x_out_shape[0]): - xi = x[i, ..., 0] - dxi = dx[i, ..., 0] - result = odl_op.derivative(xi).adjoint(dxi) - out[i] = np.asarray(result)[..., None] - - # TODO: Improve - scale = odl_op.domain.cell_volume / odl_op.range.cell_volume - out *= scale - return out - - with ops.name_scope(name + 'Grad', "ODLTensorflowLayerGrad", [x, dx]): - result = py_func(_impl, - [x, dx], - [tf.float32], - stateful=False) - - # We must manually set the output shape since tensorflow cannot - # figure it out - result = result[0] - result.set_shape(out_shape) - return result - - def tensorflow_layer_grad(op, grad): - """Thin wrapper for the gradient.""" - x = op.inputs[0] - return tensorflow_layer_grad_impl(x, grad) - - def tensorflow_layer(x): - """Implementation of the tensorflow call. - - Gradient in tensorflow is equivalent to the adjoint of the derivative - in ODL. - """ - x_shape = x.get_shape() - try: - n_x = int(x_shape[0]) - fixed_size = True - except TypeError: - n_x = x_shape[0] - fixed_size = False - - in_shape = (n_x,) + odl_op.domain.shape + (1,) - out_shape = (n_x,) + odl_op.range.shape + (1,) - - # Validate input shape - assert x_shape[1:] == odl_op.domain.shape + (1,) - - def _impl(x): - if fixed_size: - x_out_shape = out_shape - assert x.shape == in_shape - else: - x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == in_shape[1:] - - out = np.empty(x_out_shape, odl_op.range.dtype) - for i in range(x_out_shape[0]): - out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] - return out - - if differentiable: - grad = tensorflow_layer_grad - else: - grad = None - - with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: - result = py_func(_impl, - [x], - [tf.float32], - name=name_call, - stateful=False, - grad=grad) - - # We must manually set the output shape since tensorflow cannot - # figure it out - result = result[0] - result.set_shape(out_shape) - return result - - return tensorflow_layer - if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From 8817c052d5a1de18798a12b205bd08244ec8fade Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 16 Mar 2017 15:01:11 +0100 Subject: [PATCH 13/70] ENH: Add TensorflowSpace --- odl/contrib/tensorflow.py | 214 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 odl/contrib/tensorflow.py diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py new file mode 100644 index 00000000000..b544b4bb3f3 --- /dev/null +++ b/odl/contrib/tensorflow.py @@ -0,0 +1,214 @@ +import odl +import tensorflow as tf +import numpy as np + + +__all__ = ('as_tensorflow_layer', 'TensorflowSpace') + + +def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): + """Convert ``Operator`` to tensorflow layer. + + Parameters + ---------- + odl_op : `Operator` + The operator that should be wrapped to a tensorflow layer. + name : str + Tensorflow name of the operator + differentiable : boolean + True if the operator should be differentiable, otherwise assumes that + the derivative is everywhere zero. + + Returns + ------- + tensorflow_layer : callable + Callable that, when called with an `tensorflow.Tensor` of shape + `(n, *odl_op.domain.shape, 1)` returns a tensor of shape + `(n, *odl_op.range.shape, 1)` where ``n`` is the number of batches. + Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. + The `dtype` of the tensor is the same as the respective ODL spaces. + + If ``odl_op`` implements `Operator.derivative` which in turn implements + `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and + gradients propagate as expected. + """ + import tensorflow as tf + from tensorflow.python.framework import ops + + def py_func(func, inp, Tout, stateful=True, name=None, grad=None): + """Define custom py_func which takes also a grad op as argument.""" + if grad is None: + return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + else: + # Need to generate a unique name to avoid duplicates: + rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8)) + + tf.RegisterGradient(rnd_name)(grad) + g = tf.get_default_graph() + + if stateful: + override_name = 'PyFunc' + else: + override_name = 'PyFuncStateless' + + with g.gradient_override_map({override_name: rnd_name}): + return tf.py_func(func, inp, Tout, stateful=stateful, + name=name) + + def tensorflow_layer_grad_impl(x, dx): + """Implementation of the tensorflow gradient. + + Gradient in tensorflow is equivalent to the adjoint of the derivative + in ODL. + """ + x_shape = x.get_shape() + dx_shape = dx.get_shape() + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.range.shape + (1,) + out_shape = (n_x,) + odl_op.domain.shape + (1,) + + # Validate input shape + assert x_shape[1:] == odl_op.domain.shape + (1,) + assert dx_shape[1:] == odl_op.range.shape + (1,) + + def _impl(x, dx): + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape + assert dx.shape == out_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == out_shape[1:] + assert dx.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.domain.dtype) + for i in range(x_out_shape[0]): + xi = x[i, ..., 0] + dxi = dx[i, ..., 0] + result = odl_op.derivative(xi).adjoint(dxi) + out[i] = np.asarray(result)[..., None] + + # TODO: Improve + scale = odl_op.domain.cell_volume / odl_op.range.cell_volume + out *= scale + return out + + with ops.name_scope(name + 'Grad', "ODLTensorflowLayerGrad", [x, dx]): + result = py_func(_impl, + [x, dx], + [tf.float32], + stateful=False) + + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result + + def tensorflow_layer_grad(op, grad): + """Thin wrapper for the gradient.""" + x = op.inputs[0] + return tensorflow_layer_grad_impl(x, grad) + + def tensorflow_layer(x): + """Implementation of the tensorflow call. + + Gradient in tensorflow is equivalent to the adjoint of the derivative + in ODL. + """ + x_shape = x.get_shape() + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.domain.shape + (1,) + out_shape = (n_x,) + odl_op.range.shape + (1,) + + # Validate input shape + assert x_shape[1:] == odl_op.domain.shape + (1,) + + def _impl(x): + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.range.dtype) + for i in range(x_out_shape[0]): + out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] + return out + + if differentiable: + grad = tensorflow_layer_grad + else: + grad = None + + with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: + result = py_func(_impl, + [x], + [tf.float32], + name=name_call, + stateful=False, + grad=grad) + + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result + + return tensorflow_layer + + +class TensorflowSpace(odl.LinearSpace): + """A space of tensorflow Tensors.""" + def __init__(self): + odl.LinearSpace.__init__(self, odl.RealNumbers()) + + def _lincomb(self, a, x1, b, x2, out): + out.data = a * x1.data + b * x2.data + + def element(self, inp=None): + if inp in self: + return inp + elif inp is None: + return TensorflowSpaceElement(self, + tf.constant(0, dtype=tf.float32)) + else: + return TensorflowSpaceElement(self, inp) + + def __eq__(self, other): + return isinstance(other, TensorflowSpace) + + def __repr__(self): + return 'TensorflowSpace()' + + +class TensorflowSpaceElement(odl.LinearSpaceElement): + def __init__(self, space, data): + odl.LinearSpaceElement.__init__(self, space) + self.data = data + + @property + def data(self): + return self._data + + @data.setter + def data(self, value): + if isinstance(value, TensorflowSpaceElement): + raise Exception(value.data) + self._data = value + + def __repr__(self): + return '{}.element({})'.format(self.space, self.data) From 96ed70c865406ee178fc6d9b2d47860ab40b6cca Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 17 Mar 2017 09:07:02 +0100 Subject: [PATCH 14/70] BUG: Fix bug with default LinearSpace.zeros --- odl/set/space.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/set/space.py b/odl/set/space.py index d862d961975..c251deabc5a 100644 --- a/odl/set/space.py +++ b/odl/set/space.py @@ -140,7 +140,7 @@ def zero(self): """Return the zero (additive unit) element of this space.""" # Default implementation using lincomb tmp = self.element() - self.lincomb(tmp, 0, tmp, 0, tmp) + self.lincomb(0, tmp, 0, tmp, tmp) return tmp def __contains__(self, other): From 8a61714bfde8ad0fe667d92c4f328eae1a6cff8e Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 17 Mar 2017 09:09:08 +0100 Subject: [PATCH 15/70] ENH: Improvements to tensorflow support --- odl/contrib/tensorflow.py | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index b544b4bb3f3..7643c856d14 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -3,7 +3,7 @@ import numpy as np -__all__ = ('as_tensorflow_layer', 'TensorflowSpace') +__all__ = ('as_tensorflow_layer', 'TensorflowSpace', 'TensorflowSpaceOperator') def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): @@ -173,8 +173,10 @@ def _impl(x): class TensorflowSpace(odl.LinearSpace): """A space of tensorflow Tensors.""" - def __init__(self): + def __init__(self, shape): odl.LinearSpace.__init__(self, odl.RealNumbers()) + self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) else si for si in shape) + self.init_shape = tuple(si if si.value is not None else tf.Dimension(1) for si in self.shape) def _lincomb(self, a, x1, b, x2, out): out.data = a * x1.data + b * x2.data @@ -184,15 +186,20 @@ def element(self, inp=None): return inp elif inp is None: return TensorflowSpaceElement(self, - tf.constant(0, dtype=tf.float32)) + tf.zeros(self.init_shape, + dtype=tf.float32)) else: return TensorflowSpaceElement(self, inp) + def one(self): + return self.element(tf.ones(self.init_shape, + dtype=tf.float32)) + def __eq__(self, other): - return isinstance(other, TensorflowSpace) + return isinstance(other, TensorflowSpace) and other.shape == self.shape def __repr__(self): - return 'TensorflowSpace()' + return 'TensorflowSpace({})'.format(self.shape) class TensorflowSpaceElement(odl.LinearSpaceElement): @@ -212,3 +219,21 @@ def data(self, value): def __repr__(self): return '{}.element({})'.format(self.space, self.data) + + +class TensorflowSpaceOperator(odl.Operator): + def __init__(self, domain, range, func, adjoint=None, linear=False): + odl.Operator.__init__(self, domain, range, linear) + self.func = func + self.adjoint_func = adjoint + + def _call(self, x): + return self.func(x.data) + + @property + def adjoint(self): + return TensorflowSpaceOperator(self.range, + self.domain, + self.adjoint_func, + self.func, + self.is_linear) From 7977ff00de932f8f5806040e4870e8e5937bb9ac Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 17 Mar 2017 13:03:10 +0100 Subject: [PATCH 16/70] ENH: Minor performance optimization to douglas rachford solver --- odl/solvers/nonsmooth/douglas_rachford.py | 1 - 1 file changed, 1 deletion(-) diff --git a/odl/solvers/nonsmooth/douglas_rachford.py b/odl/solvers/nonsmooth/douglas_rachford.py index 01958ddc88e..81334bd07d2 100644 --- a/odl/solvers/nonsmooth/douglas_rachford.py +++ b/odl/solvers/nonsmooth/douglas_rachford.py @@ -121,7 +121,6 @@ def douglas_rachford_pd(x, f, g, L, tau, sigma, niter, composite and parallel-sum type monotone operators*. SIAM Journal on Optimization, 23.4 (2013), pp 2541--2565. """ - # Problem size m = len(L) From 2e7cc2a65d6a35b35a45bf1850628fce943ea925 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 17 Mar 2017 14:14:51 +0100 Subject: [PATCH 17/70] ENH: Add name and name scopes to tf space --- odl/contrib/tensorflow.py | 47 +++++++++++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 7643c856d14..97f43bbccf8 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -173,27 +173,60 @@ def _impl(x): class TensorflowSpace(odl.LinearSpace): """A space of tensorflow Tensors.""" - def __init__(self, shape): + def __init__(self, shape, name='ODLTensorflowSpace'): odl.LinearSpace.__init__(self, odl.RealNumbers()) self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) else si for si in shape) self.init_shape = tuple(si if si.value is not None else tf.Dimension(1) for si in self.shape) + self.name = name def _lincomb(self, a, x1, b, x2, out): - out.data = a * x1.data + b * x2.data + with tf.name_scope('{}_lincomb'.format(self.name)): + if x1 is x2: + # x1 is aligned with x2 -> out = (a+b)*x1 + out.data = (a + b) * x1.data + elif out is x1 and out is x2: + # All the vectors are aligned -> out = (a+b)*out + if (a + b) != 1: + out.data *= (a + b) + elif out is x1: + # out is aligned with x1 -> out = a*out + b*x2 + out.data = a * out.data + b * x2.data + elif out is x2: + # out is aligned with x2 -> out = a*x1 + b*out + out.data = a * x1.data + b * out.data + else: + # We have exhausted all alignment options, so x1 != x2 != out + # We now optimize for various values of a and b + if b == 0: + if a == 0: # Zero assignment -> out = 0 + out.data *= 0 + else: # Scaled copy -> out = a*x1 + out.data = a * x1.data + else: + if a == 0: # Scaled copy -> out = b*x2 + out.data = b * x2.data + elif a == 1: # No scaling in x1 -> out = x1 + b*x2 + out.data = x1.data + b * x2.data + else: # Generic case -> out = a*x1 + b*x2 + out.data = a * x1.data + b * x2.data def element(self, inp=None): if inp in self: return inp elif inp is None: - return TensorflowSpaceElement(self, - tf.zeros(self.init_shape, - dtype=tf.float32)) + return self.zero() else: return TensorflowSpaceElement(self, inp) + def zero(self): + with tf.name_scope('{}_zero'.format(self.name)): + return self.element(tf.zeros(self.init_shape, + dtype=tf.float32)) + def one(self): - return self.element(tf.ones(self.init_shape, - dtype=tf.float32)) + with tf.name_scope('{}_one'.format(self.name)): + return self.element(tf.ones(self.init_shape, + dtype=tf.float32)) def __eq__(self, other): return isinstance(other, TensorflowSpace) and other.shape == self.shape From e33285460138a4b98cf9bda992eff35f1371173c Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 13:22:29 +0100 Subject: [PATCH 18/70] ENH: Add (initial) support for tensorflow interconnectivity --- .../operator/tensorflow_cnn_tomography.py | 30 +++++++++++ .../operator/tensorflow_layer_tomography.py | 32 +++++++++++ odl/operator/oputils.py | 53 ++++++++++++++++++- 3 files changed, 114 insertions(+), 1 deletion(-) create mode 100644 examples/operator/tensorflow_cnn_tomography.py create mode 100644 examples/operator/tensorflow_layer_tomography.py diff --git a/examples/operator/tensorflow_cnn_tomography.py b/examples/operator/tensorflow_cnn_tomography.py new file mode 100644 index 00000000000..aecd8152f3b --- /dev/null +++ b/examples/operator/tensorflow_cnn_tomography.py @@ -0,0 +1,30 @@ +"""Example of how to convert an ODL operator to a tensorflow layer.""" + +import tensorflow as tf +import numpy as np +import odl + +with tf.Session() as sess: + tf.global_variables_initializer().run() + + space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space) + ray_transform = odl.tomo.RayTransform(space, geometry) + + x = tf.constant(np.asarray(ray_transform.domain.one())) + z = tf.constant(np.asarray(ray_transform.range.one())) + + odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') + y = odl_op_layer(x) + + # Evaluate using tensorflow + print(y.eval()) + + # Compare result with pure ODL + print(ray_transform(x.eval())) + + # Evaluate the adjoint of the derivative, called gradient in tensorflow + print(tf.gradients(y, [x], z)[0].eval()) + + # Compare result with pure ODL + print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/examples/operator/tensorflow_layer_tomography.py b/examples/operator/tensorflow_layer_tomography.py new file mode 100644 index 00000000000..66ad0e2124c --- /dev/null +++ b/examples/operator/tensorflow_layer_tomography.py @@ -0,0 +1,32 @@ +"""Example of how to convert an ODL ray transform to a tensorflow layer.""" + +import tensorflow as tf +import numpy as np +import odl + + +with tf.Session() as sess: + tf.global_variables_initializer().run() + + space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], + dtype='float32') + geometry = odl.tomo.parallel_beam_geometry(space) + ray_transform = odl.tomo.RayTransform(space, geometry) + + x = tf.constant(np.asarray(ray_transform.domain.one())) + z = tf.constant(np.asarray(ray_transform.range.one())) + + odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') + y = odl_op_layer(x) + + # Evaluate using tensorflow + print(y.eval()) + + # Compare result with pure ODL + print(ray_transform(x.eval())) + + # Evaluate the adjoint of the derivative, called gradient in tensorflow + print(tf.gradients(y, [x], z)[0].eval()) + + # Compare result with pure ODL + print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index a16cde45e09..e2bfc607a5b 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -21,7 +21,8 @@ from odl.util import as_flat_array __all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator', - 'as_scipy_functional', 'as_proximal_lang_operator') + 'as_scipy_functional', 'as_proximal_lang_operator', + 'as_tensorflow_layer') def matrix_representation(op): @@ -440,6 +441,56 @@ def adjoint(inp, out): norm_bound=norm_bound) +<<<<<<< 2e7cc2a65d6a35b35a45bf1850628fce943ea925 +======= +def as_tensorflow_layer(odl_op, name): + """Convert ``Operator`` to tensorflow layer.""" + import tensorflow as tf + from tensorflow.python.framework import ops + + def py_func(func, inp, Tout, stateful=True, name=None, grad=None): + """Define custom py_func which takes also a grad op as argument.""" + if grad is None: + return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + else: + # Need to generate a unique name to avoid duplicates: + rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8)) + + tf.RegisterGradient(rnd_name)(grad) + g = tf.get_default_graph() + with g.gradient_override_map({"PyFunc": rnd_name}): + return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + + def tensorflow_layer_grad_impl(x, dx): + def _impl(x, dx): + return np.asarray(odl_op.derivative(x).adjoint(dx)) + + with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): + result = py_func(_impl, + [x, dx], + [tf.float32]) + return result[0] + + # Actual gradient: + def tensorflow_layer_grad(op, grad): + x = op.inputs[0] + return tensorflow_layer_grad_impl(x, grad) + + def tensorflow_layer(x): + def _impl(x): + return np.asarray(odl_op(x)) + + with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: + result = py_func(_impl, + [x], + [tf.float32], + name=name_call, + grad=tensorflow_layer_grad) + return result[0] + + return tensorflow_layer + +>>>>>>> ENH: Add (initial) support for tensorflow interconnectivity if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From 8dba85c49d6f20c2550a88a5da7aebe2717474ea Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 13:40:30 +0100 Subject: [PATCH 19/70] ENH: Add more tensorflow examples --- .../operator/tensorflow_cnn_tomography.py | 30 ----------------- .../operator/tensorflow_layer_tomography.py | 32 ------------------- 2 files changed, 62 deletions(-) delete mode 100644 examples/operator/tensorflow_cnn_tomography.py delete mode 100644 examples/operator/tensorflow_layer_tomography.py diff --git a/examples/operator/tensorflow_cnn_tomography.py b/examples/operator/tensorflow_cnn_tomography.py deleted file mode 100644 index aecd8152f3b..00000000000 --- a/examples/operator/tensorflow_cnn_tomography.py +++ /dev/null @@ -1,30 +0,0 @@ -"""Example of how to convert an ODL operator to a tensorflow layer.""" - -import tensorflow as tf -import numpy as np -import odl - -with tf.Session() as sess: - tf.global_variables_initializer().run() - - space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space) - ray_transform = odl.tomo.RayTransform(space, geometry) - - x = tf.constant(np.asarray(ray_transform.domain.one())) - z = tf.constant(np.asarray(ray_transform.range.one())) - - odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') - y = odl_op_layer(x) - - # Evaluate using tensorflow - print(y.eval()) - - # Compare result with pure ODL - print(ray_transform(x.eval())) - - # Evaluate the adjoint of the derivative, called gradient in tensorflow - print(tf.gradients(y, [x], z)[0].eval()) - - # Compare result with pure ODL - print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/examples/operator/tensorflow_layer_tomography.py b/examples/operator/tensorflow_layer_tomography.py deleted file mode 100644 index 66ad0e2124c..00000000000 --- a/examples/operator/tensorflow_layer_tomography.py +++ /dev/null @@ -1,32 +0,0 @@ -"""Example of how to convert an ODL ray transform to a tensorflow layer.""" - -import tensorflow as tf -import numpy as np -import odl - - -with tf.Session() as sess: - tf.global_variables_initializer().run() - - space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], - dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space) - ray_transform = odl.tomo.RayTransform(space, geometry) - - x = tf.constant(np.asarray(ray_transform.domain.one())) - z = tf.constant(np.asarray(ray_transform.range.one())) - - odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') - y = odl_op_layer(x) - - # Evaluate using tensorflow - print(y.eval()) - - # Compare result with pure ODL - print(ray_transform(x.eval())) - - # Evaluate the adjoint of the derivative, called gradient in tensorflow - print(tf.gradients(y, [x], z)[0].eval()) - - # Compare result with pure ODL - print(ray_transform.derivative(x.eval()).adjoint(z.eval())) From 2586f50184e4a2723d398a386375697f5ecb8a44 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 17:27:01 +0100 Subject: [PATCH 20/70] ENH: Improve handling of shapes of tensorflow layers --- odl/operator/oputils.py | 86 +++++++++++++++++++++++++++++++++++------ 1 file changed, 75 insertions(+), 11 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index e2bfc607a5b..284b6ff0069 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -441,10 +441,29 @@ def adjoint(inp, out): norm_bound=norm_bound) -<<<<<<< 2e7cc2a65d6a35b35a45bf1850628fce943ea925 -======= -def as_tensorflow_layer(odl_op, name): - """Convert ``Operator`` to tensorflow layer.""" +def as_tensorflow_layer(odl_op, name='ODLOperator'): + """Convert ``Operator`` to tensorflow layer. + + Parameters + ---------- + odl_op : `Operator` + The operator that should be wrapped to a tensorflow layer. + name : str + Tensorflow name of the operator + + Returns + ------- + tensorflow_layer : callable + Callable that, when called with an `tensorflow.Tensor` of shape + `(n, *odl_op.domain.shape, 1)` returns a tensor of shape + `(n, *odl_op.range.shape, 1)` where ``n`` is the number of batches. + Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. + The `dtype` of the tensor is the same as the respective ODL spaces. + + If ``odl_op`` implements `Operator.derivative` which in turn implements + `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and + gradients propagate as expected. + """ import tensorflow as tf from tensorflow.python.framework import ops @@ -459,26 +478,66 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): tf.RegisterGradient(rnd_name)(grad) g = tf.get_default_graph() with g.gradient_override_map({"PyFunc": rnd_name}): - return tf.py_func(func, inp, Tout, stateful=stateful, name=name) + return tf.py_func(func, inp, Tout, stateful=stateful, + name=name) def tensorflow_layer_grad_impl(x, dx): + """Implementation of the tensorflow gradient. + + Gradient in tensorflow is equivalent to the adjoint of the derivative + in ODL. + """ + x_shape = x.get_shape() + dx_shape = dx.get_shape() + out_shape = (int(x_shape[0]),) + odl_op.domain.shape + (1,) + + # Validate input shape + assert x_shape[0] == dx_shape[0] + assert x_shape[1:] == odl_op.domain.shape + (1,) + assert dx_shape[1:] == odl_op.range.shape + (1,) + def _impl(x, dx): - return np.asarray(odl_op.derivative(x).adjoint(dx)) + out = np.empty(out_shape, odl_op.domain.dtype) + for i in range(x_shape[0]): + xi = x[i, ..., 0] + dxi = dx[i, ..., 0] + result = odl_op.derivative(xi).adjoint(dxi) + out[i] = np.asarray(result)[..., None] + return out with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): result = py_func(_impl, [x, dx], [tf.float32]) - return result[0] - # Actual gradient: + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result + def tensorflow_layer_grad(op, grad): + """Thin wrapper for the gradient.""" x = op.inputs[0] return tensorflow_layer_grad_impl(x, grad) def tensorflow_layer(x): + """Implementation of the tensorflow call. + + Gradient in tensorflow is equivalent to the adjoint of the derivative + in ODL. + """ + x_shape = x.get_shape() + out_shape = (int(x_shape[0]),) + odl_op.range.shape + (1,) + + # Validate input shape + assert x_shape[1:] == odl_op.domain.shape + (1,) + def _impl(x): - return np.asarray(odl_op(x)) + out = np.empty(out_shape, odl_op.range.dtype) + for i in range(x_shape[0]): + out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] + return out with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: result = py_func(_impl, @@ -486,11 +545,16 @@ def _impl(x): [tf.float32], name=name_call, grad=tensorflow_layer_grad) - return result[0] + + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result return tensorflow_layer ->>>>>>> ENH: Add (initial) support for tensorflow interconnectivity + if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From 70fe54f715f7e2f8cab374215c3dd40e759128ce Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 17:27:11 +0100 Subject: [PATCH 21/70] ENH: Add more tensorflow examples --- .../operator/tensorflow_tomography_cnn_gradient.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index dcc314f1e2b..dc50e2d1a24 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -26,20 +26,17 @@ def var(x): return tf.Variable(tf.constant(x, dtype='float32')) -def create_variable(name, shape, stddev=0.01): - variable = tf.Variable(tf.truncated_normal(shape, stddev=stddev), name=name) - return variable - - with tf.Session() as sess: - # Create ODL data structures - size = 128 + size = 256 space = odl.uniform_discr([-64, -64], [64, 64], [size, size], dtype='float32') geometry = odl.tomo.parallel_beam_geometry(space) ray_transform = odl.tomo.RayTransform(space, geometry) - + phantom = odl.phantom.shepp_logan(space, True) + data = ray_transform(phantom) + noisy_data = data + odl.phantom.white_noise(ray_transform.range) * np.mean(data) * 0.1 + fbp = odl.tomo.fbp_op(ray_transform)(noisy_data) # Create tensorflow layer from odl operator odl_op_layer = odl.as_tensorflow_layer(ray_transform, @@ -128,3 +125,4 @@ def create_variable(name, shape, stddev=0.01): print('iter={}, training loss={}, validation loss={}'.format(i, loss_training, loss_value)) callback(space.element(x_value[0])) + From 7c7943b87e35db638ac27f36a40fdb9c42abd684 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 21:39:50 +0100 Subject: [PATCH 22/70] ENH: Optimizations to the tensorflow interface --- odl/operator/oputils.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 284b6ff0069..3a4a4d1002f 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -503,12 +503,17 @@ def _impl(x, dx): dxi = dx[i, ..., 0] result = odl_op.derivative(xi).adjoint(dxi) out[i] = np.asarray(result)[..., None] + + # TODO: Improve + scale = odl_op.domain.cell_volume / odl_op.range.cell_volume + out *= scale return out with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): result = py_func(_impl, [x, dx], - [tf.float32]) + [tf.float32], + stateful=False) # We must manually set the output shape since tensorflow cannot # figure it out @@ -544,6 +549,7 @@ def _impl(x): [x], [tf.float32], name=name_call, + stateful=False, grad=tensorflow_layer_grad) # We must manually set the output shape since tensorflow cannot From ccc8220dc436d4c82f4bc4d7e40ea28858b82e97 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Feb 2017 21:40:14 +0100 Subject: [PATCH 23/70] ENH: Improved CNN gradient tomography --- examples/operator/tensorflow_tomography_cnn_gradient.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index dc50e2d1a24..cd294e2ad5f 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -28,15 +28,12 @@ def var(x): with tf.Session() as sess: # Create ODL data structures - size = 256 + size = 128 space = odl.uniform_discr([-64, -64], [64, 64], [size, size], dtype='float32') geometry = odl.tomo.parallel_beam_geometry(space) ray_transform = odl.tomo.RayTransform(space, geometry) - phantom = odl.phantom.shepp_logan(space, True) - data = ray_transform(phantom) - noisy_data = data + odl.phantom.white_noise(ray_transform.range) * np.mean(data) * 0.1 - fbp = odl.tomo.fbp_op(ray_transform)(noisy_data) + # Create tensorflow layer from odl operator odl_op_layer = odl.as_tensorflow_layer(ray_transform, @@ -125,4 +122,3 @@ def var(x): print('iter={}, training loss={}, validation loss={}'.format(i, loss_training, loss_value)) callback(space.element(x_value[0])) - From 921fe33b8eb76b4623a03899aaec96c2211702e9 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Feb 2017 13:49:56 +0100 Subject: [PATCH 24/70] ENH: Support variable size input for tensorflow --- odl/operator/oputils.py | 42 +++++++++++++++++++++++++++++++++++------ 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 3a4a4d1002f..47be47010de 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -489,7 +489,15 @@ def tensorflow_layer_grad_impl(x, dx): """ x_shape = x.get_shape() dx_shape = dx.get_shape() - out_shape = (int(x_shape[0]),) + odl_op.domain.shape + (1,) + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.range.shape + (1,) + out_shape = (n_x,) + odl_op.domain.shape + (1,) # Validate input shape assert x_shape[0] == dx_shape[0] @@ -497,8 +505,15 @@ def tensorflow_layer_grad_impl(x, dx): assert dx_shape[1:] == odl_op.range.shape + (1,) def _impl(x, dx): - out = np.empty(out_shape, odl_op.domain.dtype) - for i in range(x_shape[0]): + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.domain.dtype) + for i in range(x_out_shape[0]): xi = x[i, ..., 0] dxi = dx[i, ..., 0] result = odl_op.derivative(xi).adjoint(dxi) @@ -533,14 +548,29 @@ def tensorflow_layer(x): in ODL. """ x_shape = x.get_shape() - out_shape = (int(x_shape[0]),) + odl_op.range.shape + (1,) + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.domain.shape + (1,) + out_shape = (n_x,) + odl_op.range.shape + (1,) # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) def _impl(x): - out = np.empty(out_shape, odl_op.range.dtype) - for i in range(x_shape[0]): + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.range.dtype) + for i in range(x_out_shape[0]): out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] return out From ce186b5253d7533f25374d36779347b49bd915fc Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Feb 2017 14:27:05 +0100 Subject: [PATCH 25/70] ENH: Improved cnn gradient example --- .../operator/tensorflow_tomography_cnn_gradient.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index cd294e2ad5f..f9d29368c5d 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -26,7 +26,13 @@ def var(x): return tf.Variable(tf.constant(x, dtype='float32')) +def create_variable(name, shape): + variable = tf.Variable(tf.truncated_normal(shape, stddev=0.01), name=name) + return variable + + with tf.Session() as sess: + # Create ODL data structures size = 128 space = odl.uniform_discr([-64, -64], [64, 64], [size, size], @@ -67,13 +73,13 @@ def var(x): # Create constant right hand side - w1 = tf.Variable(tf.truncated_normal([3, 3, 7, 32], stddev=0.01)) + w1 = create_variable('w1', shape=[3, 3, 7, 32]) b1 = var(np.ones(32) * 0.1) - w2 = tf.Variable(tf.truncated_normal([3, 3, 32, 32], stddev=0.01)) + w2 = create_variable('w2', shape=[3, 3, 32, 32]) b2 = var(np.ones(32) * 0.1) - w3 = tf.Variable(tf.truncated_normal([3, 3, 32, 6], stddev=0.01)) + w3 = create_variable('w3', shape=[3, 3, 32, 6]) b3 = var(np.random.randn(6) * 0.001) n_iter = 5 @@ -108,7 +114,7 @@ def var(x): # Solve with an ODL callback to see what happens callback = odl.solvers.CallbackShow() - for i in range(10000): + for i in range(1000): _, loss_training = sess.run([optimizer, loss], feed_dict={learning_rate: 0.001, x_0: x_arr[1:], From bf896c22750e451f9ecfd57f22d0d06b931de204 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Feb 2017 15:38:56 +0100 Subject: [PATCH 26/70] ENH: Improvements to the CNN tensorflow example --- .../operator/tensorflow_tomography_cnn_gradient.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/operator/tensorflow_tomography_cnn_gradient.py index f9d29368c5d..dcc314f1e2b 100644 --- a/examples/operator/tensorflow_tomography_cnn_gradient.py +++ b/examples/operator/tensorflow_tomography_cnn_gradient.py @@ -26,8 +26,8 @@ def var(x): return tf.Variable(tf.constant(x, dtype='float32')) -def create_variable(name, shape): - variable = tf.Variable(tf.truncated_normal(shape, stddev=0.01), name=name) +def create_variable(name, shape, stddev=0.01): + variable = tf.Variable(tf.truncated_normal(shape, stddev=stddev), name=name) return variable @@ -73,13 +73,13 @@ def create_variable(name, shape): # Create constant right hand side - w1 = create_variable('w1', shape=[3, 3, 7, 32]) + w1 = tf.Variable(tf.truncated_normal([3, 3, 7, 32], stddev=0.01)) b1 = var(np.ones(32) * 0.1) - w2 = create_variable('w2', shape=[3, 3, 32, 32]) + w2 = tf.Variable(tf.truncated_normal([3, 3, 32, 32], stddev=0.01)) b2 = var(np.ones(32) * 0.1) - w3 = create_variable('w3', shape=[3, 3, 32, 6]) + w3 = tf.Variable(tf.truncated_normal([3, 3, 32, 6], stddev=0.01)) b3 = var(np.random.randn(6) * 0.001) n_iter = 5 @@ -114,7 +114,7 @@ def create_variable(name, shape): # Solve with an ODL callback to see what happens callback = odl.solvers.CallbackShow() - for i in range(1000): + for i in range(10000): _, loss_training = sess.run([optimizer, loss], feed_dict={learning_rate: 0.001, x_0: x_arr[1:], From 4a34da3d82ab9c8a75fd8adfdb4bc3fb21e5c045 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 3 Mar 2017 17:12:27 +0100 Subject: [PATCH 27/70] BUG: Fix bug with derivative of tensorflow operators --- odl/operator/oputils.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 47be47010de..782fd8de14d 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -477,7 +477,13 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): tf.RegisterGradient(rnd_name)(grad) g = tf.get_default_graph() - with g.gradient_override_map({"PyFunc": rnd_name}): + + if stateful: + override_name = 'PyFunc' + else: + override_name = 'PyFuncStateless' + + with g.gradient_override_map({override_name: rnd_name}): return tf.py_func(func, inp, Tout, stateful=stateful, name=name) @@ -500,7 +506,6 @@ def tensorflow_layer_grad_impl(x, dx): out_shape = (n_x,) + odl_op.domain.shape + (1,) # Validate input shape - assert x_shape[0] == dx_shape[0] assert x_shape[1:] == odl_op.domain.shape + (1,) assert dx_shape[1:] == odl_op.range.shape + (1,) @@ -508,9 +513,11 @@ def _impl(x, dx): if fixed_size: x_out_shape = out_shape assert x.shape == in_shape + assert dx.shape == out_shape else: x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == in_shape[1:] + assert x.shape[1:] == out_shape[1:] + assert dx.shape[1:] == in_shape[1:] out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): @@ -524,7 +531,7 @@ def _impl(x, dx): out *= scale return out - with ops.name_scope(name, "ODLTensorflowLayerGrad", [x, dx]): + with ops.name_scope(name + 'Grad', "ODLTensorflowLayerGrad", [x, dx]): result = py_func(_impl, [x, dx], [tf.float32], From b2414531c589a988c473ffc5cb7c9bebcb2b631c Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 13 Mar 2017 17:41:42 +0100 Subject: [PATCH 28/70] ENH: Allow non-differentialbe operators in tensorflow --- odl/operator/oputils.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 782fd8de14d..b928cbf8de3 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -441,7 +441,7 @@ def adjoint(inp, out): norm_bound=norm_bound) -def as_tensorflow_layer(odl_op, name='ODLOperator'): +def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): """Convert ``Operator`` to tensorflow layer. Parameters @@ -450,6 +450,9 @@ def as_tensorflow_layer(odl_op, name='ODLOperator'): The operator that should be wrapped to a tensorflow layer. name : str Tensorflow name of the operator + differentiable : boolean + True if the operator should be differentiable, otherwise assumes that + the derivative is everywhere zero. Returns ------- @@ -581,13 +584,18 @@ def _impl(x): out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] return out + if differentiable: + grad = tensorflow_layer_grad + else: + grad = None + with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: result = py_func(_impl, [x], [tf.float32], name=name_call, stateful=False, - grad=tensorflow_layer_grad) + grad=grad) # We must manually set the output shape since tensorflow cannot # figure it out From 21fb1e1126a0df9c8d4f9c43786d3b0a9c267c73 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 16 Mar 2017 15:00:57 +0100 Subject: [PATCH 29/70] ENH: Add odl.contrib for less-than-optimal code --- odl/operator/oputils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index b928cbf8de3..1582c3ce8b7 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -21,8 +21,7 @@ from odl.util import as_flat_array __all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator', - 'as_scipy_functional', 'as_proximal_lang_operator', - 'as_tensorflow_layer') + 'as_scipy_functional', 'as_proximal_lang_operator') def matrix_representation(op): @@ -441,6 +440,7 @@ def adjoint(inp, out): norm_bound=norm_bound) +<<<<<<< b2414531c589a988c473ffc5cb7c9bebcb2b631c def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): """Convert ``Operator`` to tensorflow layer. From 4241f0b9c6beef4e8895a2325214592ce5b56c17 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 16 Mar 2017 15:01:11 +0100 Subject: [PATCH 30/70] ENH: Add TensorflowSpace --- odl/contrib/tensorflow.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 97f43bbccf8..171a014eebc 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -173,6 +173,7 @@ def _impl(x): class TensorflowSpace(odl.LinearSpace): """A space of tensorflow Tensors.""" + def __init__(self, shape, name='ODLTensorflowSpace'): odl.LinearSpace.__init__(self, odl.RealNumbers()) self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) else si for si in shape) @@ -222,7 +223,7 @@ def zero(self): with tf.name_scope('{}_zero'.format(self.name)): return self.element(tf.zeros(self.init_shape, dtype=tf.float32)) - + def one(self): with tf.name_scope('{}_one'.format(self.name)): return self.element(tf.ones(self.init_shape, From 148a8ecf0badf451e3452b993f3996901282ac6a Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 17 Mar 2017 14:14:51 +0100 Subject: [PATCH 31/70] ENH: Add name and name scopes to tf space --- odl/contrib/tensorflow.py | 1 + 1 file changed, 1 insertion(+) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 171a014eebc..d8865132e3d 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -172,6 +172,7 @@ def _impl(x): class TensorflowSpace(odl.LinearSpace): + """A space of tensorflow Tensors.""" def __init__(self, shape, name='ODLTensorflowSpace'): From 0c3f2d39e032c9c574ad97129a3433d475f6f65e Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 20 Mar 2017 18:23:51 +0100 Subject: [PATCH 32/70] tmp --- odl/contrib/tensorflow.py | 55 ++++++++++++++++++++++++++++++++------- 1 file changed, 45 insertions(+), 10 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index d8865132e3d..39bad5389de 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -70,14 +70,21 @@ def tensorflow_layer_grad_impl(x, dx): n_x = x_shape[0] fixed_size = False - in_shape = (n_x,) + odl_op.range.shape + (1,) + if odl_op.is_functional: + in_shape = (n_x, 1) + else: + in_shape = (n_x,) + odl_op.range.shape + (1,) out_shape = (n_x,) + odl_op.domain.shape + (1,) # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) - assert dx_shape[1:] == odl_op.range.shape + (1,) + if odl_op.is_functional: + assert dx_shape[1:] == (1,) + else: + assert dx_shape[1:] == odl_op.range.shape + (1,) def _impl(x, dx): + if fixed_size: x_out_shape = out_shape assert x.shape == in_shape @@ -89,13 +96,28 @@ def _impl(x, dx): out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): - xi = x[i, ..., 0] - dxi = dx[i, ..., 0] - result = odl_op.derivative(xi).adjoint(dxi) - out[i] = np.asarray(result)[..., None] + if odl_op.is_functional: + xi = x[i, ..., 0] + dxi = dx[i, 0] + out[i, ..., 0] = np.asarray(odl_op.gradient(xi)) * dxi + else: + xi = x[i, ..., 0] + dxi = dx[i, ..., 0] + result = odl_op.derivative(xi).adjoint(dxi) + out[i, ..., 0] = np.asarray(result) # TODO: Improve - scale = odl_op.domain.cell_volume / odl_op.range.cell_volume + try: + dom_weight = odl_op.domain.weighting.const + except AttributeError: + dom_weight = 1.0 + + try: + ran_weight = odl_op.range.weighting.const + except AttributeError: + ran_weight = 1.0 + + scale = dom_weight / ran_weight out *= scale return out @@ -131,22 +153,35 @@ def tensorflow_layer(x): fixed_size = False in_shape = (n_x,) + odl_op.domain.shape + (1,) - out_shape = (n_x,) + odl_op.range.shape + (1,) + if odl_op.is_functional: + out_shape = (n_x, 1) + else: + out_shape = (n_x,) + odl_op.range.shape + (1,) + # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) def _impl(x): + print('CALL ASSERT') if fixed_size: x_out_shape = out_shape assert x.shape == in_shape else: x_out_shape = (x.shape[0],) + out_shape[1:] assert x.shape[1:] == in_shape[1:] + print('CALL ASSERT OK') - out = np.empty(x_out_shape, odl_op.range.dtype) + out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): - out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] + print('CALL LOOP') + if odl_op.is_functional: + out[i, 0] = odl_op(x[i, ..., 0]) + else: + out[i, ..., 0] = np.asarray(odl_op(x[i, ..., 0])) + print('CALL DONE') + print(name) + print('shape {}'.format(out.shape)) return out if differentiable: From cb6790146f150f7812ef2e5768ff2c80d846c1d9 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 20 Mar 2017 18:48:03 +0100 Subject: [PATCH 33/70] ENH: Add derivative of gradient of l1 norm --- odl/solvers/functional/default_functionals.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/odl/solvers/functional/default_functionals.py b/odl/solvers/functional/default_functionals.py index 51cfdaa94d8..8e1ee525bd7 100644 --- a/odl/solvers/functional/default_functionals.py +++ b/odl/solvers/functional/default_functionals.py @@ -144,6 +144,10 @@ def _call(self, x): """Apply the gradient operator to the given point.""" return x.ufuncs.sign() + def derivative(self, x): + """Derivative is a.e. zero.""" + return ZeroOperator(self.domain) + return L1Gradient() elif self.exponent == 2: class L2Gradient(Operator): From 5a78b97ae2538e419baf368118330be490b36b10 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 20 Mar 2017 18:48:17 +0100 Subject: [PATCH 34/70] ENH: Enable functionals with tensorflow --- odl/contrib/tensorflow.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 39bad5389de..318171177ed 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -84,7 +84,6 @@ def tensorflow_layer_grad_impl(x, dx): assert dx_shape[1:] == odl_op.range.shape + (1,) def _impl(x, dx): - if fixed_size: x_out_shape = out_shape assert x.shape == in_shape @@ -158,30 +157,24 @@ def tensorflow_layer(x): else: out_shape = (n_x,) + odl_op.range.shape + (1,) - # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) def _impl(x): - print('CALL ASSERT') if fixed_size: x_out_shape = out_shape assert x.shape == in_shape else: x_out_shape = (x.shape[0],) + out_shape[1:] assert x.shape[1:] == in_shape[1:] - print('CALL ASSERT OK') out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): - print('CALL LOOP') if odl_op.is_functional: out[i, 0] = odl_op(x[i, ..., 0]) else: out[i, ..., 0] = np.asarray(odl_op(x[i, ..., 0])) - print('CALL DONE') - print(name) - print('shape {}'.format(out.shape)) + return out if differentiable: From 03c3f34a48554eb4aab65c5c7bbaf170d8f0aee0 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 21 Mar 2017 09:13:37 +0100 Subject: [PATCH 35/70] BUG: Fix bug with tensorflow layer gradient with fixed sizes --- odl/contrib/tensorflow.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 318171177ed..1474ad962aa 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -86,8 +86,8 @@ def tensorflow_layer_grad_impl(x, dx): def _impl(x, dx): if fixed_size: x_out_shape = out_shape - assert x.shape == in_shape - assert dx.shape == out_shape + assert x.shape == out_shape + assert dx.shape == in_shape else: x_out_shape = (x.shape[0],) + out_shape[1:] assert x.shape[1:] == out_shape[1:] From ac977a32f75d5e86ba1c2a67328579ed4fb87b87 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 21 Mar 2017 12:13:48 +0100 Subject: [PATCH 36/70] API: Change shape of tensors returned from functionals with tensorflow --- odl/contrib/tensorflow.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 1474ad962aa..f1d4441581a 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -71,7 +71,7 @@ def tensorflow_layer_grad_impl(x, dx): fixed_size = False if odl_op.is_functional: - in_shape = (n_x, 1) + in_shape = (n_x,) else: in_shape = (n_x,) + odl_op.range.shape + (1,) out_shape = (n_x,) + odl_op.domain.shape + (1,) @@ -79,7 +79,7 @@ def tensorflow_layer_grad_impl(x, dx): # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) if odl_op.is_functional: - assert dx_shape[1:] == (1,) + assert dx_shape[1:] == () else: assert dx_shape[1:] == odl_op.range.shape + (1,) @@ -97,7 +97,7 @@ def _impl(x, dx): for i in range(x_out_shape[0]): if odl_op.is_functional: xi = x[i, ..., 0] - dxi = dx[i, 0] + dxi = dx[i] out[i, ..., 0] = np.asarray(odl_op.gradient(xi)) * dxi else: xi = x[i, ..., 0] @@ -153,7 +153,7 @@ def tensorflow_layer(x): in_shape = (n_x,) + odl_op.domain.shape + (1,) if odl_op.is_functional: - out_shape = (n_x, 1) + out_shape = (n_x,) else: out_shape = (n_x,) + odl_op.range.shape + (1,) @@ -171,7 +171,7 @@ def _impl(x): out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): if odl_op.is_functional: - out[i, 0] = odl_op(x[i, ..., 0]) + out[i] = odl_op(x[i, ..., 0]) else: out[i, ..., 0] = np.asarray(odl_op(x[i, ..., 0])) From 762aa11fcb1705b56a34dbb6bca6b766390a6e08 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 23 Mar 2017 12:27:24 +0100 Subject: [PATCH 37/70] ENH: Improve doc of tensorflow layer --- odl/contrib/tensorflow.py | 84 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 81 insertions(+), 3 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index f1d4441581a..d7a288d61c5 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -14,7 +14,7 @@ def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): odl_op : `Operator` The operator that should be wrapped to a tensorflow layer. name : str - Tensorflow name of the operator + Tensorflow name of the operator. differentiable : boolean True if the operator should be differentiable, otherwise assumes that the derivative is everywhere zero. @@ -84,6 +84,43 @@ def tensorflow_layer_grad_impl(x, dx): assert dx_shape[1:] == odl_op.range.shape + (1,) def _impl(x, dx): + """Implementation of the adjoint of the derivative. + + Returns :: + + odl_op.derivative(x).adjoint(dx) + + Parameters + ---------- + x : `numpy.ndarray` + Point(s) in which the derivative should be taken. + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. This is a constant if ``fixed_size`` is true, + otherwise it is dynamic. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. + dx : `numpy.ndarray` + Point(s) in which the adjoint of the derivative of the operator + should be evaluated. + The axes are: + 0 : batch id. Should be pairwise matched with ``x``. + 1, ..., M-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + + Returns + ------- + result : `numpy.ndarray` + Result of the computation. + + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. + """ if fixed_size: x_out_shape = out_shape assert x.shape == out_shape @@ -140,8 +177,28 @@ def tensorflow_layer_grad(op, grad): def tensorflow_layer(x): """Implementation of the tensorflow call. - Gradient in tensorflow is equivalent to the adjoint of the derivative - in ODL. + Returns a `tensorflow.Tensor` that represents a lazy application of + ``odl_op`` to ``x``. + + Parameters + ---------- + x : `tensorflow.Tensor` + Point(s) to which the layer should be applied. + The axes are: + 0 : batch id. Can be fixed or dynamic. + 1, ..., M-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + + Returns + ------- + result : `tensorflow.Tensor` + Lazy result of the computation. + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. """ x_shape = x.get_shape() try: @@ -161,6 +218,27 @@ def tensorflow_layer(x): assert x_shape[1:] == odl_op.domain.shape + (1,) def _impl(x): + """Implementation of the tensorflow layer. + + Parameters + ---------- + x : `numpy.ndarray` + Point(s) in which the operator should be evaluated. + The axes are: + 0 : batch id. This is a constant if ``fixed_size`` is true, + otherwise it is dynamic. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + + Returns + ------- + result : `numpy.ndarray` + Result of the computation. + The axes are: + 0 : batch id. Data is pairwise matched with ``x``. + 1, ..., M-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + """ if fixed_size: x_out_shape = out_shape assert x.shape == in_shape From 6809e118b546bc434d52f8105edf482274a315a1 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 24 Mar 2017 10:58:07 +0100 Subject: [PATCH 38/70] ENH: Fix bug with astra cuda in case of several callers --- odl/tomo/backends/astra_cuda.py | 121 +++++++++++++++++--------------- 1 file changed, 66 insertions(+), 55 deletions(-) diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 6869d77b8eb..9666f1defa3 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -29,6 +29,7 @@ from odl.tomo.geometry import ( Geometry, Parallel2dGeometry, FanFlatGeometry, Parallel3dAxisGeometry, ConeFlatGeometry) +from multiprocessing import Lock __all__ = ('ASTRA_CUDA_AVAILABLE', @@ -62,6 +63,10 @@ def __init__(self, geometry, reco_space, proj_space): self.create_ids() + # Create a mutually exclusive lock so that two callers cant use the + # same shared resource at the same time. + self.mutex = Lock() + def call_forward(self, vol_data, out=None): """Run an ASTRA forward projection on the given data using the GPU. @@ -79,36 +84,37 @@ def call_forward(self, vol_data, out=None): Projection data resulting from the application of the projector. If ``out`` was provided, the returned object is a reference to it. """ - assert vol_data in self.reco_space - if out is not None: - assert out in self.proj_space - else: - out = self.proj_space.element() - - # Copy data to GPU memory - if self.geometry.ndim == 2: - astra.data2d.store(self.vol_id, vol_data.asarray()) - elif self.geometry.ndim == 3: - astra.data3d.store(self.vol_id, vol_data.asarray()) - else: - raise RuntimeError('unknown ndim') - - # Run algorithm - astra.algorithm.run(self.algo_id) - - # Copy result to host - if self.geometry.ndim == 2: - out[:] = self.out_array - elif self.geometry.ndim == 3: - out[:] = np.swapaxes(self.out_array, 0, 1).reshape( - self.proj_space.shape) - - # Fix scaling to weight by pixel size - if isinstance(self.geometry, Parallel2dGeometry): - # parallel2d scales with pixel stride - out *= 1 / float(self.geometry.det_partition.cell_volume) - - return out + with self.mutex: + assert vol_data in self.reco_space + if out is not None: + assert out in self.proj_space + else: + out = self.proj_space.element() + + # Copy data to GPU memory + if self.geometry.ndim == 2: + astra.data2d.store(self.vol_id, vol_data.asarray()) + elif self.geometry.ndim == 3: + astra.data3d.store(self.vol_id, vol_data.asarray()) + else: + raise RuntimeError('unknown ndim') + + # Run algorithm + astra.algorithm.run(self.algo_id) + + # Copy result to host + if self.geometry.ndim == 2: + out[:] = self.out_array + elif self.geometry.ndim == 3: + out[:] = np.swapaxes(self.out_array, 0, 1).reshape( + self.proj_space.shape) + + # Fix scaling to weight by pixel size + if isinstance(self.geometry, Parallel2dGeometry): + # parallel2d scales with pixel stride + out *= 1 / float(self.geometry.det_partition.cell_volume) + + return out def create_ids(self): """Create ASTRA objects.""" @@ -205,6 +211,10 @@ def __init__(self, geometry, reco_space, proj_space): self.proj_space = proj_space self.create_ids() + # Create a mutually exclusive lock so that two callers cant use the + # same shared resource at the same time. + self.mutex = Lock() + def call_backward(self, proj_data, out=None): """Run an ASTRA back-projection on the given data using the GPU. @@ -223,33 +233,34 @@ def call_backward(self, proj_data, out=None): back-projector. If ``out`` was provided, the returned object is a reference to it. """ - assert proj_data in self.proj_space - if out is not None: - assert out in self.reco_space - else: - out = self.reco_space.element() - - # Copy data to GPU memory - if self.geometry.ndim == 2: - astra.data2d.store(self.sino_id, proj_data.asarray()) - elif self.geometry.ndim == 3: - shape = (-1,) + self.geometry.det_partition.shape - reshaped_proj_data = proj_data.asarray().reshape(shape) - swapped_proj_data = np.ascontiguousarray( - np.swapaxes(reshaped_proj_data, 0, 1)) - astra.data3d.store(self.sino_id, swapped_proj_data) - - # Run algorithm - astra.algorithm.run(self.algo_id) - - # Copy result to CPU memory - out[:] = self.out_array + with self.mutex: + assert proj_data in self.proj_space + if out is not None: + assert out in self.reco_space + else: + out = self.reco_space.element() + + # Copy data to GPU memory + if self.geometry.ndim == 2: + astra.data2d.store(self.sino_id, proj_data.asarray()) + elif self.geometry.ndim == 3: + shape = (-1,) + self.geometry.det_partition.shape + reshaped_proj_data = proj_data.asarray().reshape(shape) + swapped_proj_data = np.ascontiguousarray( + np.swapaxes(reshaped_proj_data, 0, 1)) + astra.data3d.store(self.sino_id, swapped_proj_data) + + # Run algorithm + astra.algorithm.run(self.algo_id) + + # Copy result to CPU memory + out[:] = self.out_array - # Fix scaling to weight by pixel/voxel size - out *= astra_cuda_bp_scaling_factor( - self.proj_space, self.reco_space, self.geometry) + # Fix scaling to weight by pixel/voxel size + out *= astra_cuda_bp_scaling_factor( + self.proj_space, self.reco_space, self.geometry) - return out + return out def create_ids(self): """Create ASTRA objects.""" From fdd41dd3f4d442b7c7ed328cb97b3d1176b2d06a Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 24 Mar 2017 11:18:28 +0100 Subject: [PATCH 39/70] ENH: Improve doc and add name scopes to tensorflow --- odl/contrib/tensorflow.py | 421 +++++++++++++++++++++----------------- 1 file changed, 229 insertions(+), 192 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index d7a288d61c5..816c3ccc48f 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -6,15 +6,16 @@ __all__ = ('as_tensorflow_layer', 'TensorflowSpace', 'TensorflowSpaceOperator') -def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): +def as_tensorflow_layer(odl_op, default_name='ODLOperator', + differentiable=True): """Convert ``Operator`` to tensorflow layer. Parameters ---------- odl_op : `Operator` The operator that should be wrapped to a tensorflow layer. - name : str - Tensorflow name of the operator. + default_name : str + Default name for tensorflow layers created. differentiable : boolean True if the operator should be differentiable, otherwise assumes that the derivative is everywhere zero. @@ -40,141 +41,166 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): if grad is None: return tf.py_func(func, inp, Tout, stateful=stateful, name=name) else: - # Need to generate a unique name to avoid duplicates: - rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8)) - - tf.RegisterGradient(rnd_name)(grad) - g = tf.get_default_graph() - if stateful: override_name = 'PyFunc' else: override_name = 'PyFuncStateless' + # Need to generate a unique name to avoid duplicates: + rnd_name = override_name + 'Grad' + str(np.random.randint(0, 1E+8)) + + tf.RegisterGradient(rnd_name)(grad) + g = tf.get_default_graph() + with g.gradient_override_map({override_name: rnd_name}): return tf.py_func(func, inp, Tout, stateful=stateful, name=name) - def tensorflow_layer_grad_impl(x, dx): + def tensorflow_layer_grad_impl(x, dx, name): """Implementation of the tensorflow gradient. Gradient in tensorflow is equivalent to the adjoint of the derivative in ODL. - """ - x_shape = x.get_shape() - dx_shape = dx.get_shape() - try: - n_x = int(x_shape[0]) - fixed_size = True - except TypeError: - n_x = x_shape[0] - fixed_size = False - - if odl_op.is_functional: - in_shape = (n_x,) - else: - in_shape = (n_x,) + odl_op.range.shape + (1,) - out_shape = (n_x,) + odl_op.domain.shape + (1,) - # Validate input shape - assert x_shape[1:] == odl_op.domain.shape + (1,) - if odl_op.is_functional: - assert dx_shape[1:] == () - else: - assert dx_shape[1:] == odl_op.range.shape + (1,) - - def _impl(x, dx): - """Implementation of the adjoint of the derivative. - - Returns :: - - odl_op.derivative(x).adjoint(dx) - - Parameters - ---------- - x : `numpy.ndarray` - Point(s) in which the derivative should be taken. - If ``odl_op`` is an `Operator` the axes are: - 0 : batch id. This is a constant if ``fixed_size`` is true, - otherwise it is dynamic. - 1, ..., N-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. - If ``odl_op`` is a `Functional` the axes are: - 0 : batch id. - dx : `numpy.ndarray` - Point(s) in which the adjoint of the derivative of the operator - should be evaluated. - The axes are: - 0 : batch id. Should be pairwise matched with ``x``. - 1, ..., M-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. - - Returns - ------- - result : `numpy.ndarray` - Result of the computation. - - If ``odl_op`` is an `Operator` the axes are: - 0 : batch id. - 1, ..., N-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. - If ``odl_op`` is a `Functional` the axes are: - 0 : batch id. - """ - if fixed_size: - x_out_shape = out_shape - assert x.shape == out_shape - assert dx.shape == in_shape - else: - x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == out_shape[1:] - assert dx.shape[1:] == in_shape[1:] - - out = np.empty(x_out_shape, odl_op.domain.dtype) - for i in range(x_out_shape[0]): - if odl_op.is_functional: - xi = x[i, ..., 0] - dxi = dx[i] - out[i, ..., 0] = np.asarray(odl_op.gradient(xi)) * dxi - else: - xi = x[i, ..., 0] - dxi = dx[i, ..., 0] - result = odl_op.derivative(xi).adjoint(dxi) - out[i, ..., 0] = np.asarray(result) + Returns a `tensorflow.Tensor` that represents a lazy application of :: - # TODO: Improve - try: - dom_weight = odl_op.domain.weighting.const - except AttributeError: - dom_weight = 1.0 + odl_op.derivative(x).adjoint(dx) + + Parameters + ---------- + x : `tensorflow.Tensor` + Point(s) to which the layer should be applied. + The axes are: + 0 : batch id. Can be fixed or dynamic. + 1, ..., M-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + name : string + Name of the tensor. + Returns + ------- + result : `tensorflow.Tensor` + Lazy result of the computation. + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. + """ + with tf.name_scope(name): + x_shape = x.get_shape() + dx_shape = dx.get_shape() try: - ran_weight = odl_op.range.weighting.const - except AttributeError: - ran_weight = 1.0 - - scale = dom_weight / ran_weight - out *= scale - return out - - with ops.name_scope(name + 'Grad', "ODLTensorflowLayerGrad", [x, dx]): - result = py_func(_impl, - [x, dx], - [tf.float32], - stateful=False) - - # We must manually set the output shape since tensorflow cannot - # figure it out - result = result[0] - result.set_shape(out_shape) - return result - - def tensorflow_layer_grad(op, grad): - """Thin wrapper for the gradient.""" - x = op.inputs[0] - return tensorflow_layer_grad_impl(x, grad) - - def tensorflow_layer(x): + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + if odl_op.is_functional: + in_shape = (n_x,) + else: + in_shape = (n_x,) + odl_op.range.shape + (1,) + out_shape = (n_x,) + odl_op.domain.shape + (1,) + + # Validate input shape + assert x_shape[1:] == odl_op.domain.shape + (1,) + if odl_op.is_functional: + assert dx_shape[1:] == () + else: + assert dx_shape[1:] == odl_op.range.shape + (1,) + + def _impl(x, dx): + """Implementation of the adjoint of the derivative. + + Returns :: + + odl_op.derivative(x).adjoint(dx) + + Parameters + ---------- + x : `numpy.ndarray` + Point(s) in which the derivative should be taken. + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. This is a constant if ``fixed_size`` is + true, otherwise it is dynamic. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. + dx : `numpy.ndarray` + Point(s) in which the adjoint of the derivative of the + operator should be evaluated. + The axes are: + 0 : batch id. Should be pairwise matched with ``x``. + 1, ..., M-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + + Returns + ------- + result : `numpy.ndarray` + Result of the computation. + + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. + """ + if fixed_size: + x_out_shape = out_shape + assert x.shape == out_shape + assert dx.shape == in_shape + else: + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == out_shape[1:] + assert dx.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.domain.dtype) + out[:] = np.nan + for i in range(x_out_shape[0]): + if odl_op.is_functional: + xi = x[i, ..., 0] + dxi = dx[i] + out[i, ..., 0] = np.asarray(odl_op.gradient(xi)) * dxi + else: + xi = x[i, ..., 0] + dxi = dx[i, ..., 0] + result = odl_op.derivative(xi).adjoint(dxi) + out[i, ..., 0] = np.asarray(result) + + # TODO: Improve + try: + dom_weight = odl_op.domain.weighting.const + except AttributeError: + dom_weight = 1.0 + + try: + ran_weight = odl_op.range.weighting.const + except AttributeError: + ran_weight = 1.0 + + scale = dom_weight / ran_weight + out *= scale + + return out + + with ops.name_scope(name + '_pyfunc', values=[x, dx]) as name_call: + result = py_func(_impl, + [x, dx], + [tf.float32], + name=name_call, + stateful=True) + + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result + + def tensorflow_layer(x, name=None): """Implementation of the tensorflow call. Returns a `tensorflow.Tensor` that represents a lazy application of @@ -188,6 +214,8 @@ def tensorflow_layer(x): 0 : batch id. Can be fixed or dynamic. 1, ..., M-2 : spatial dimensions of data. n-1 : (currently) unused data channel. + name : string + Name of the tensor. Default: Defaultname. Returns ------- @@ -200,79 +228,88 @@ def tensorflow_layer(x): If ``odl_op`` is a `Functional` the axes are: 0 : batch id. """ - x_shape = x.get_shape() - try: - n_x = int(x_shape[0]) - fixed_size = True - except TypeError: - n_x = x_shape[0] - fixed_size = False - - in_shape = (n_x,) + odl_op.domain.shape + (1,) - if odl_op.is_functional: - out_shape = (n_x,) - else: - out_shape = (n_x,) + odl_op.range.shape + (1,) - - # Validate input shape - assert x_shape[1:] == odl_op.domain.shape + (1,) - - def _impl(x): - """Implementation of the tensorflow layer. - - Parameters - ---------- - x : `numpy.ndarray` - Point(s) in which the operator should be evaluated. - The axes are: - 0 : batch id. This is a constant if ``fixed_size`` is true, - otherwise it is dynamic. - 1, ..., N-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. - - Returns - ------- - result : `numpy.ndarray` - Result of the computation. - The axes are: - 0 : batch id. Data is pairwise matched with ``x``. - 1, ..., M-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. - """ - if fixed_size: - x_out_shape = out_shape - assert x.shape == in_shape - else: - x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == in_shape[1:] + if name is None: + name = default_name - out = np.empty(x_out_shape, odl_op.domain.dtype) - for i in range(x_out_shape[0]): - if odl_op.is_functional: - out[i] = odl_op(x[i, ..., 0]) + with tf.name_scope(name): + x_shape = x.get_shape() + try: + n_x = int(x_shape[0]) + fixed_size = True + except TypeError: + n_x = x_shape[0] + fixed_size = False + + in_shape = (n_x,) + odl_op.domain.shape + (1,) + if odl_op.is_functional: + out_shape = (n_x,) + else: + out_shape = (n_x,) + odl_op.range.shape + (1,) + + # Validate input shape + assert x_shape[1:] == odl_op.domain.shape + (1,) + + def _impl(x): + """Implementation of the tensorflow layer. + + Parameters + ---------- + x : `numpy.ndarray` + Point(s) in which the operator should be evaluated. + The axes are: + 0 : batch id. This is a constant if ``fixed_size`` is + true, otherwise it is dynamic. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + + Returns + ------- + result : `numpy.ndarray` + Result of the computation. + The axes are: + 0 : batch id. Data is pairwise matched with ``x``. + 1, ..., M-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + """ + if fixed_size: + x_out_shape = out_shape + assert x.shape == in_shape else: - out[i, ..., 0] = np.asarray(odl_op(x[i, ..., 0])) - - return out - - if differentiable: - grad = tensorflow_layer_grad - else: - grad = None - - with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: - result = py_func(_impl, - [x], - [tf.float32], - name=name_call, - stateful=False, - grad=grad) - - # We must manually set the output shape since tensorflow cannot - # figure it out - result = result[0] - result.set_shape(out_shape) - return result + x_out_shape = (x.shape[0],) + out_shape[1:] + assert x.shape[1:] == in_shape[1:] + + out = np.empty(x_out_shape, odl_op.domain.dtype) + out[:] = np.nan + for i in range(x_out_shape[0]): + if odl_op.is_functional: + out[i] = odl_op(x[i, ..., 0]) + else: + out[i, ..., 0] = np.asarray(odl_op(x[i, ..., 0])) + + return out + + if differentiable: + def tensorflow_layer_grad(op, grad): + """Thin wrapper for the gradient.""" + x = op.inputs[0] + return tensorflow_layer_grad_impl(x, grad, + name=name + '_grad') + else: + tensorflow_layer_grad = None + + with ops.name_scope(name + '_pyfunc', values=[x]) as name_call: + result = py_func(_impl, + [x], + [tf.float32], + name=name_call, + stateful=True, + grad=tensorflow_layer_grad) + + # We must manually set the output shape since tensorflow cannot + # figure it out + result = result[0] + result.set_shape(out_shape) + return result return tensorflow_layer From daf7d83464d092bd9672491b39017d0f38abea86 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 24 Mar 2017 13:12:48 +0100 Subject: [PATCH 40/70] MAINT: Bug fixes and test improvements --- odl/contrib/tensorflow.py | 8 ++++++-- odl/tomo/backends/astra_cuda.py | 6 +++--- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 816c3ccc48f..78f04b9db71 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -320,8 +320,12 @@ class TensorflowSpace(odl.LinearSpace): def __init__(self, shape, name='ODLTensorflowSpace'): odl.LinearSpace.__init__(self, odl.RealNumbers()) - self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) else si for si in shape) - self.init_shape = tuple(si if si.value is not None else tf.Dimension(1) for si in self.shape) + self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) + else si + for si in shape) + self.init_shape = tuple(si if si.value is not None + else tf.Dimension(1) + for si in self.shape) self.name = name def _lincomb(self, a, x1, b, x2, out): diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 9666f1defa3..63d09d5fab8 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -211,9 +211,9 @@ def __init__(self, geometry, reco_space, proj_space): self.proj_space = proj_space self.create_ids() - # Create a mutually exclusive lock so that two callers cant use the - # same shared resource at the same time. - self.mutex = Lock() + # Create a mutually exclusive lock so that two callers cant use the + # same shared resource at the same time. + self.mutex = Lock() def call_backward(self, proj_data, out=None): """Run an ASTRA back-projection on the given data using the GPU. From 889430d9f65bdd40e5d505f0c93e3e947a89b852 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 24 Mar 2017 16:11:34 +0100 Subject: [PATCH 41/70] ENH: Change ellipse_phantom to use radians --- odl/phantom/geometric.py | 10 ++++++---- odl/phantom/transmission.py | 10 ++++++---- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/odl/phantom/geometric.py b/odl/phantom/geometric.py index 7c7b9da25b3..b40a675dc98 100644 --- a/odl/phantom/geometric.py +++ b/odl/phantom/geometric.py @@ -353,7 +353,7 @@ def _ellipse_phantom_2d(space, ellipses): b_squared = ellip[2] ** 2 x0 = ellip[3] y0 = ellip[4] - theta = ellip[5] * np.pi / 180 + theta = ellip[5] scales = [1 / a_squared, 1 / b_squared] center = (np.array([x0, y0]) + 1.0) / 2.0 @@ -474,9 +474,9 @@ def _ellipsoid_phantom_3d(space, ellipsoids): x0 = ellip[4] y0 = ellip[5] z0 = ellip[6] - phi = ellip[7] * np.pi / 180 - theta = ellip[8] * np.pi / 180 - psi = ellip[9] * np.pi / 180 + phi = ellip[7] + theta = ellip[8] + psi = ellip[9] scales = [1 / a_squared, 1 / b_squared, 1 / c_squared] center = (np.array([x0, y0, z0]) + 1.0) / 2.0 @@ -563,6 +563,8 @@ def ellipsoid_phantom(space, ellipsoids): The ellipsoids need to be given such that the ellipsoids fall in the rectangle [-1, -1] x [1, 1] or equivalent in 3d. + The angles are given in radians. + Notes ----- The phantom is created by adding the values of each ellipse. The ellipses diff --git a/odl/phantom/transmission.py b/odl/phantom/transmission.py index f605accf840..e3c3a4bf04f 100644 --- a/odl/phantom/transmission.py +++ b/odl/phantom/transmission.py @@ -27,11 +27,12 @@ def _shepp_logan_ellipse_2d(): This assumes that the ellipses are contained in the square [-1, -1]x[-1, -1]. """ + rad18 = np.deg2rad(18.0) # value axisx axisy x y rotation return [[2.00, .6900, .9200, 0.0000, 0.0000, 0], [-.98, .6624, .8740, 0.0000, -.0184, 0], - [-.02, .1100, .3100, 0.2200, 0.0000, -18], - [-.02, .1600, .4100, -.2200, 0.0000, 18], + [-.02, .1100, .3100, 0.2200, 0.0000, -rad18], + [-.02, .1600, .4100, -.2200, 0.0000, rad18], [0.01, .2100, .2500, 0.0000, 0.3500, 0], [0.01, .0460, .0460, 0.0000, 0.1000, 0], [0.01, .0460, .0460, 0.0000, -.1000, 0], @@ -46,11 +47,12 @@ def _shepp_logan_ellipsoids_3d(): This assumes that the ellipses are contained in the cube [-1, -1, -1]x[1, 1, 1]. """ + rad18 = np.deg2rad(18.0) # value axisx axisy axisz, x y z rotation return [[2.00, .6900, .9200, .810, 0.0000, 0.0000, 0.00, 0.0, 0, 0], [-.98, .6624, .8740, .780, 0.0000, -.0184, 0.00, 0.0, 0, 0], - [-.02, .1100, .3100, .220, 0.2200, 0.0000, 0.00, -18, 0, 0], - [-.02, .1600, .4100, .280, -.2200, 0.0000, 0.00, 18., 0, 0], + [-.02, .1100, .3100, .220, 0.2200, 0.0000, 0.00, -rad18, 0, 0], + [-.02, .1600, .4100, .280, -.2200, 0.0000, 0.00, rad18, 0, 0], [0.01, .2100, .2500, .410, 0.0000, 0.3500, 0.00, 0.0, 0, 0], [0.01, .0460, .0460, .050, 0.0000, 0.1000, 0.00, 0.0, 0, 0], [0.01, .0460, .0460, .050, 0.0000, -.1000, 0.00, 0.0, 0, 0], From 0a09fb82f685538b58b0eb0c048792b4215bccd0 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 27 Mar 2017 13:24:09 +0200 Subject: [PATCH 42/70] ENH: Fix return dtypes in tensorflow space --- odl/contrib/tensorflow.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 78f04b9db71..6ddf9e2545f 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -190,7 +190,7 @@ def _impl(x, dx): with ops.name_scope(name + '_pyfunc', values=[x, dx]) as name_call: result = py_func(_impl, [x, dx], - [tf.float32], + [odl_op.domain.dtype], name=name_call, stateful=True) @@ -249,6 +249,9 @@ def tensorflow_layer(x, name=None): # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) + out_dtype = getattr(odl_op.range, 'dtype', + odl_op.domain.dtype) + def _impl(x): """Implementation of the tensorflow layer. @@ -278,7 +281,7 @@ def _impl(x): x_out_shape = (x.shape[0],) + out_shape[1:] assert x.shape[1:] == in_shape[1:] - out = np.empty(x_out_shape, odl_op.domain.dtype) + out = np.empty(x_out_shape, out_dtype) out[:] = np.nan for i in range(x_out_shape[0]): if odl_op.is_functional: @@ -300,7 +303,7 @@ def tensorflow_layer_grad(op, grad): with ops.name_scope(name + '_pyfunc', values=[x]) as name_call: result = py_func(_impl, [x], - [tf.float32], + [out_dtype], name=name_call, stateful=True, grad=tensorflow_layer_grad) From 616d9142611b6c4a9c93b13802ea08a4e3c9b32f Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Wed, 29 Mar 2017 11:03:50 +0200 Subject: [PATCH 43/70] ENH: Remove unncessary warning when calling matplotlib.pause --- odl/util/graphics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/odl/util/graphics.py b/odl/util/graphics.py index 0037aa04c75..17476884ae2 100644 --- a/odl/util/graphics.py +++ b/odl/util/graphics.py @@ -32,7 +32,6 @@ def warning_free_pause(): message="Using default event loop until " "function specific to this GUI is " "implemented") - plt.pause(0.0001) From a619f4c3c57c27d7c57e4ee17deddded71ce465b Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Sun, 2 Apr 2017 16:55:40 +0200 Subject: [PATCH 44/70] MAINT: Move tensorflow examples to contrib folder --- .../{operator => contrib/tensorflow}/tensorflow_layer_matrix.py | 0 .../tensorflow}/tensorflow_layer_ray_transform.py | 0 .../{operator => contrib/tensorflow}/tensorflow_tomography.py | 0 .../{operator => contrib/tensorflow}/tensorflow_tomography_cnn.py | 0 .../tensorflow}/tensorflow_tomography_cnn_gradient.py | 0 5 files changed, 0 insertions(+), 0 deletions(-) rename examples/{operator => contrib/tensorflow}/tensorflow_layer_matrix.py (100%) rename examples/{operator => contrib/tensorflow}/tensorflow_layer_ray_transform.py (100%) rename examples/{operator => contrib/tensorflow}/tensorflow_tomography.py (100%) rename examples/{operator => contrib/tensorflow}/tensorflow_tomography_cnn.py (100%) rename examples/{operator => contrib/tensorflow}/tensorflow_tomography_cnn_gradient.py (100%) diff --git a/examples/operator/tensorflow_layer_matrix.py b/examples/contrib/tensorflow/tensorflow_layer_matrix.py similarity index 100% rename from examples/operator/tensorflow_layer_matrix.py rename to examples/contrib/tensorflow/tensorflow_layer_matrix.py diff --git a/examples/operator/tensorflow_layer_ray_transform.py b/examples/contrib/tensorflow/tensorflow_layer_ray_transform.py similarity index 100% rename from examples/operator/tensorflow_layer_ray_transform.py rename to examples/contrib/tensorflow/tensorflow_layer_ray_transform.py diff --git a/examples/operator/tensorflow_tomography.py b/examples/contrib/tensorflow/tensorflow_tomography.py similarity index 100% rename from examples/operator/tensorflow_tomography.py rename to examples/contrib/tensorflow/tensorflow_tomography.py diff --git a/examples/operator/tensorflow_tomography_cnn.py b/examples/contrib/tensorflow/tensorflow_tomography_cnn.py similarity index 100% rename from examples/operator/tensorflow_tomography_cnn.py rename to examples/contrib/tensorflow/tensorflow_tomography_cnn.py diff --git a/examples/operator/tensorflow_tomography_cnn_gradient.py b/examples/contrib/tensorflow/tensorflow_tomography_cnn_gradient.py similarity index 100% rename from examples/operator/tensorflow_tomography_cnn_gradient.py rename to examples/contrib/tensorflow/tensorflow_tomography_cnn_gradient.py From 5baaee6eb21d99f490ba74a54a1a60169af4ff4e Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Sun, 2 Apr 2017 17:36:07 +0200 Subject: [PATCH 45/70] ENH: Improve doc of tensorflow support --- odl/contrib/tensorflow.py | 102 ++++++++++++++++++++++++++++++++------ 1 file changed, 87 insertions(+), 15 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 6ddf9e2545f..894cd19c601 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -1,6 +1,37 @@ +# Copyright 2014-2017 The ODL development group +# +# This file is part of ODL. +# +# ODL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ODL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with ODL. If not, see . + +"""Utilities for interfacing ODL with Tensorflow.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from builtins import range, str +from future import standard_library +standard_library.install_aliases() + import odl -import tensorflow as tf import numpy as np +import uuid +try: + # TODO: Lazy import to improve performance + import tensorflow as tf + TENSORFLOW_AVAILABLE = True +except: + TENSORFLOW_AVAILABLE = False __all__ = ('as_tensorflow_layer', 'TensorflowSpace', 'TensorflowSpaceOperator') @@ -17,7 +48,7 @@ def as_tensorflow_layer(odl_op, default_name='ODLOperator', default_name : str Default name for tensorflow layers created. differentiable : boolean - True if the operator should be differentiable, otherwise assumes that + True if the layer should be differentiable, otherwise assumes that the derivative is everywhere zero. Returns @@ -33,11 +64,33 @@ def as_tensorflow_layer(odl_op, default_name='ODLOperator', `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and gradients propagate as expected. """ - import tensorflow as tf + if not TENSORFLOW_AVAILABLE: + raise ImportError('Tensorflow not installed') + from tensorflow.python.framework import ops def py_func(func, inp, Tout, stateful=True, name=None, grad=None): - """Define custom py_func which takes also a grad op as argument.""" + """Define custom py_func which takes also a grad op as argument. + + We need to overwrite this function since the default tensorflow py_func + does not support custom gradients. + + Parameters + ---------- + func : callable + Python function that takes and returns numpy arrays. + inp : sequence of `tensorflow.Tensor` + Input tensors for the function + Tout : sequence of `tensorflow.dtype` + Datatype of the output(s). + stateful : bool + If the function has internal state, i.e. if calling the function + with a given input always gives the same output. + name : string + Name of the python function + grad : callbable + Gradient of the function. + """ if grad is None: return tf.py_func(func, inp, Tout, stateful=stateful, name=name) else: @@ -47,7 +100,7 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): override_name = 'PyFuncStateless' # Need to generate a unique name to avoid duplicates: - rnd_name = override_name + 'Grad' + str(np.random.randint(0, 1E+8)) + rnd_name = override_name + 'Grad' + str(uuid.uuid4()) tf.RegisterGradient(rnd_name)(grad) g = tf.get_default_graph() @@ -68,10 +121,20 @@ def tensorflow_layer_grad_impl(x, dx, name): Parameters ---------- - x : `tensorflow.Tensor` - Point(s) to which the layer should be applied. + x : `numpy.ndarray` + Point(s) in which the derivative should be taken. + If ``odl_op`` is an `Operator` the axes are: + 0 : batch id. This is a constant if ``fixed_size`` is + true, otherwise it is dynamic. + 1, ..., N-2 : spatial dimensions of data. + n-1 : (currently) unused data channel. + If ``odl_op`` is a `Functional` the axes are: + 0 : batch id. + dx : `tensorflow.Tensor` + Point(s) in which the adjoint of the derivative of the + operator should be evaluated. The axes are: - 0 : batch id. Can be fixed or dynamic. + 0 : batch id. Should be pairwise matched with ``x``. 1, ..., M-2 : spatial dimensions of data. n-1 : (currently) unused data channel. name : string @@ -89,6 +152,7 @@ def tensorflow_layer_grad_impl(x, dx, name): 0 : batch id. """ with tf.name_scope(name): + # Validate the input/output shape x_shape = x.get_shape() dx_shape = dx.get_shape() try: @@ -104,7 +168,6 @@ def tensorflow_layer_grad_impl(x, dx, name): in_shape = (n_x,) + odl_op.range.shape + (1,) out_shape = (n_x,) + odl_op.domain.shape + (1,) - # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) if odl_op.is_functional: assert dx_shape[1:] == () @@ -149,6 +212,7 @@ def _impl(x, dx): If ``odl_op`` is a `Functional` the axes are: 0 : batch id. """ + # Validate the shape of the given input if fixed_size: x_out_shape = out_shape assert x.shape == out_shape @@ -158,8 +222,8 @@ def _impl(x, dx): assert x.shape[1:] == out_shape[1:] assert dx.shape[1:] == in_shape[1:] + # Evaluate the operator on all inputs in the batch. out = np.empty(x_out_shape, odl_op.domain.dtype) - out[:] = np.nan for i in range(x_out_shape[0]): if odl_op.is_functional: xi = x[i, ..., 0] @@ -171,7 +235,8 @@ def _impl(x, dx): result = odl_op.derivative(xi).adjoint(dxi) out[i, ..., 0] = np.asarray(result) - # TODO: Improve + # Rescale the domain/range according to the weighting since + # tensorflow does not have weighted spaces. try: dom_weight = odl_op.domain.weighting.const except AttributeError: @@ -192,7 +257,7 @@ def _impl(x, dx): [x, dx], [odl_op.domain.dtype], name=name_call, - stateful=True) + stateful=False) # We must manually set the output shape since tensorflow cannot # figure it out @@ -232,6 +297,7 @@ def tensorflow_layer(x, name=None): name = default_name with tf.name_scope(name): + # Validate input shape x_shape = x.get_shape() try: n_x = int(x_shape[0]) @@ -246,7 +312,6 @@ def tensorflow_layer(x, name=None): else: out_shape = (n_x,) + odl_op.range.shape + (1,) - # Validate input shape assert x_shape[1:] == odl_op.domain.shape + (1,) out_dtype = getattr(odl_op.range, 'dtype', @@ -274,6 +339,7 @@ def _impl(x): 1, ..., M-2 : spatial dimensions of data. n-1 : (currently) unused data channel. """ + # Validate input shape if fixed_size: x_out_shape = out_shape assert x.shape == in_shape @@ -281,8 +347,8 @@ def _impl(x): x_out_shape = (x.shape[0],) + out_shape[1:] assert x.shape[1:] == in_shape[1:] + # Evaluate the operator on all inputs in the batch. out = np.empty(x_out_shape, out_dtype) - out[:] = np.nan for i in range(x_out_shape[0]): if odl_op.is_functional: out[i] = odl_op(x[i, ..., 0]) @@ -305,7 +371,7 @@ def tensorflow_layer_grad(op, grad): [x], [out_dtype], name=name_call, - stateful=True, + stateful=False, grad=tensorflow_layer_grad) # We must manually set the output shape since tensorflow cannot @@ -388,6 +454,9 @@ def __repr__(self): class TensorflowSpaceElement(odl.LinearSpaceElement): + + """Elements in TensorflowSpace.""" + def __init__(self, space, data): odl.LinearSpaceElement.__init__(self, space) self.data = data @@ -407,6 +476,9 @@ def __repr__(self): class TensorflowSpaceOperator(odl.Operator): + + """Wrap ODL operator so that it acts on TensorflowSpace elements.""" + def __init__(self, domain, range, func, adjoint=None, linear=False): odl.Operator.__init__(self, domain, range, linear) self.func = func From b2ab020553920b018f7a699bed822745ac76e348 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Sun, 2 Apr 2017 17:48:50 +0200 Subject: [PATCH 46/70] ENH: Improve doc of tensorflow --- doc/source/guide/overview_functionality.rst | 103 ++++++++++++++++++++ odl/contrib/tensorflow.py | 8 ++ 2 files changed, 111 insertions(+) create mode 100644 doc/source/guide/overview_functionality.rst diff --git a/doc/source/guide/overview_functionality.rst b/doc/source/guide/overview_functionality.rst new file mode 100644 index 00000000000..25071a0aa0f --- /dev/null +++ b/doc/source/guide/overview_functionality.rst @@ -0,0 +1,103 @@ +.. _overview_functionality: + +######################### +Overview of functionality +######################### + + +This section gives a short overview of the functionality available in ODL. + + +What is vectorization? +====================== + +In general, :term:`vectorization` means that a function can be evaluated on a whole array of values +at once instead of looping over individual entries. This is very important for performance in an +interpreted language like python, since loops are usually very slow compared to compiled languages. + +Technically, vectorization in Numpy works through the `Universal functions (ufunc)`_ interface. It +is fast because all loops over data are implemented in C, and the resulting implementations are +exposed to Python for each function individually. + + +How to use Numpy's ufuncs? +========================== + +The easiest way to write fast custom mathematical functions in Python is to use the +`available ufuncs`_ and compose them to a new function:: + + def gaussian(x): + # Negation, powers and scaling are vectorized, of course. + return np.exp(-x ** 2 / 2) + + def step(x): + # np.where checks the condition in the first argument and + # returns the second for `True`, otherwise the third. The + # last two arguments can be arrays, too. + # Note that also the comparison operation is vectorized. + return np.where(x[0] <= 0, 0, 1) + +This should cover a very large range of useful functions already (basic arithmetic is vectorized, +too!). An even larger list of `special functions`_ are available in the Scipy package. + + +Usage in ODL +============ + +Python functions are in most cases used as input to a discretization process. For example, we may +want to discretize a two-dimensional Gaussian function: + + >>> def gaussian2(x): + ... return np.exp(-(x[0] ** 2 + x[1] ** 2) / 2) + +on the rectangle [-5, 5] x [-5, 5] with 100 pixels in each +dimension. The code for this is simply + + >>> # Note that the minimum and maxiumum coordinates are given as + >>> # vectors, not one interval at a time. + >>> discr = odl.uniform_discr([-5, -5], [5, 5], (100, 100)) + + >>> # This creates an element in the discretized space ``discr`` + >>> gaussian_discr = discr.element(gaussian2) + +What happens behind the scenes is that ``discr`` creates a :term:`discretization` object which +has a built-in method ``element`` to turn continuous functions into discrete arrays by evaluating +them at a set of grid points. In the example above, this grid is a uniform sampling of the rectangle +by 100 points per dimension. + +To make this process fast, ``element`` assumes that the function is written in a way that not only +supports vectorization, but also guarantees that the output has the correct shape. The function +receives a :term:`meshgrid` tuple as input, in the above case consisting of two vectors:: + + >>> mesh = discr.meshgrid + >>> mesh[0].shape + (100, 1) + >>> mesh[1].shape + (1, 100) + +When inserted into the function, the final shape of the output is determined by Numpy's +`broadcasting rules`_. For the Gaussian function, Numpy will conclude that the output shape must +be ``(100, 100)`` since the arrays in ``mesh`` are added after squaring. This size is the same +as expected by the discretization. + +If a function does not use all components of the input, ODL tries to broadcast the result to the shape of the discretized space:: + + >>> def gaussian_const_x0(x): + ... return np.exp(-x[1] ** 2 / 2) # no x[0] -> broadcasting + + >>> gaussian_const_x0(mesh).shape + (1, 100) + >>> discr.element(gaussian_const_x0).shape + (100, 100) + + +Further reading +=============== + +`Scipy Lecture notes on Numpy `_ + + +.. _Universal functions (ufunc): http://docs.scipy.org/doc/numpy/reference/ufuncs.html +.. _available ufuncs: http://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs +.. _special functions: http://docs.scipy.org/doc/scipy/reference/special.html +.. _broadcasting rules: http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 894cd19c601..7a0a2bda2f0 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -156,6 +156,7 @@ def tensorflow_layer_grad_impl(x, dx, name): x_shape = x.get_shape() dx_shape = dx.get_shape() try: + # Lazy check if the first dimension is dynamic n_x = int(x_shape[0]) fixed_size = True except TypeError: @@ -300,6 +301,7 @@ def tensorflow_layer(x, name=None): # Validate input shape x_shape = x.get_shape() try: + # Lazy check if the first dimension is dynamic n_x = int(x_shape[0]) fixed_size = True except TypeError: @@ -494,3 +496,9 @@ def adjoint(self): self.adjoint_func, self.func, self.is_linear) + + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() From 283031afedfd589af356901b89878f3da8eda2b8 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Sun, 2 Apr 2017 18:09:16 +0200 Subject: [PATCH 47/70] ENH: Further enhancements to tensorflow doc --- odl/contrib/tensorflow.py | 40 +++++++++++++++++++++++---------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow.py index 7a0a2bda2f0..addac92bed8 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow.py @@ -23,7 +23,9 @@ from future import standard_library standard_library.install_aliases() -import odl +from odl.set import LinearSpace, LinearSpaceElement, RealNumbers +from odl.operator import Operator + import numpy as np import uuid try: @@ -39,30 +41,36 @@ def as_tensorflow_layer(odl_op, default_name='ODLOperator', differentiable=True): - """Convert ``Operator`` to tensorflow layer. + """Convert `Operator` or `Functional` into tensorflow layer. Parameters ---------- - odl_op : `Operator` + odl_op : `Operator` or `Functional` The operator that should be wrapped to a tensorflow layer. default_name : str Default name for tensorflow layers created. differentiable : boolean - True if the layer should be differentiable, otherwise assumes that - the derivative is everywhere zero. + True if the layer should be differentiable, in which case ``odl_op`` + should implement `Operator.derivative` which in turn implements + `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and + gradients propagate as expected. + + Otherwise assumes that the derivative is everywhere zero. Returns ------- tensorflow_layer : callable Callable that, when called with an `tensorflow.Tensor` of shape - `(n, *odl_op.domain.shape, 1)` returns a tensor of shape - `(n, *odl_op.range.shape, 1)` where ``n`` is the number of batches. + `(n, *odl_op.domain.shape, 1)` where ``n`` is the number of batches + returns another `tensorflow.Tensor`. + + If ``odl_op`` is an `Operator`, the shape of the returned tensor is + `(n, *odl_op.range.shape, 1)`. + Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. + The `dtype` of the tensor is the same as the respective ODL spaces. - If ``odl_op`` implements `Operator.derivative` which in turn implements - `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and - gradients propagate as expected. """ if not TENSORFLOW_AVAILABLE: raise ImportError('Tensorflow not installed') @@ -385,12 +393,12 @@ def tensorflow_layer_grad(op, grad): return tensorflow_layer -class TensorflowSpace(odl.LinearSpace): +class TensorflowSpace(LinearSpace): """A space of tensorflow Tensors.""" def __init__(self, shape, name='ODLTensorflowSpace'): - odl.LinearSpace.__init__(self, odl.RealNumbers()) + LinearSpace.__init__(self, RealNumbers()) self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) else si for si in shape) @@ -455,12 +463,12 @@ def __repr__(self): return 'TensorflowSpace({})'.format(self.shape) -class TensorflowSpaceElement(odl.LinearSpaceElement): +class TensorflowSpaceElement(LinearSpaceElement): """Elements in TensorflowSpace.""" def __init__(self, space, data): - odl.LinearSpaceElement.__init__(self, space) + LinearSpaceElement.__init__(self, space) self.data = data @property @@ -477,12 +485,12 @@ def __repr__(self): return '{}.element({})'.format(self.space, self.data) -class TensorflowSpaceOperator(odl.Operator): +class TensorflowSpaceOperator(Operator): """Wrap ODL operator so that it acts on TensorflowSpace elements.""" def __init__(self, domain, range, func, adjoint=None, linear=False): - odl.Operator.__init__(self, domain, range, linear) + Operator.__init__(self, domain, range, linear) self.func = func self.adjoint_func = adjoint From ef163ca3d7ebfef88852a3a338603c44b1d6bc97 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 3 Apr 2017 10:11:14 +0200 Subject: [PATCH 48/70] ENH: Add util/statistics with MSE and PSNR calculators --- odl/util/__init__.py | 3 ++ odl/util/statistics.py | 103 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 106 insertions(+) create mode 100644 odl/util/statistics.py diff --git a/odl/util/__init__.py b/odl/util/__init__.py index 9f57e6dcd00..54fff755bfd 100644 --- a/odl/util/__init__.py +++ b/odl/util/__init__.py @@ -27,6 +27,9 @@ from .numerics import * __all__ += numerics.__all__ +from .statistics import * +__all__ += statistics.__all__ + from .vectorization import * __all__ += vectorization.__all__ diff --git a/odl/util/statistics.py b/odl/util/statistics.py new file mode 100644 index 00000000000..f90df7e7009 --- /dev/null +++ b/odl/util/statistics.py @@ -0,0 +1,103 @@ +# Copyright 2014-2016 The ODL development group +# +# This file is part of ODL. +# +# ODL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ODL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with ODL. If not, see . + +"""Utilities for computing statistics on images.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import numpy as np + + +__all__ = ('psnr',) + + +def mse(true, noisy): + """Return the Mean Squared Errror. + + Parameters + ---------- + true : array-like + noisy : array-like + + Returns + ------- + psnr : float + + Examples + -------- + >>> true = [1, 1, 1, 1, 1] + >>> noisy = [1, 1, 1, 1, 2] + >>> result = mse(true, noisy) + >>> print('{:.3f}'.format(result)) + 0.200 + """ + true = np.asarray(true) + noisy = np.asarray(noisy) + return np.mean((true - noisy) ** 2) + + +def psnr(true, noisy): + """Return the Peak signal-to-noise ratio. + + Parameters + ---------- + true : array-like + noisy : array-like + + Returns + ------- + psnr : float + + Examples + -------- + >>> true = [1, 1, 1, 1, 1] + >>> noisy = [1, 1, 1, 1, 2] + >>> result = psnr(true, noisy) + >>> print('{:.3f}'.format(result)) + 6.990 + + If true == noisy, the result is infinite + + >>> psnr([1, 1], [1, 1]) + inf + + If `true == 0` but `noisy != 0`, the result is negative infinity + + >>> psnr(0, 1) + -inf + """ + true = np.asarray(true) + noisy = np.asarray(noisy) + + mse_result = mse(true, noisy) + max_true = np.max(np.abs(true)) + + if mse_result == 0: + return np.inf + elif max_true == 0: + return -np.inf + else: + return 20 * np.log10(max_true) - 10 * np.log10(mse_result) + + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() From fb336c710f480146534794be551f94baa557a38b Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 13 Apr 2017 18:44:09 +0200 Subject: [PATCH 49/70] MAINT: Doc fixes to MSE --- odl/util/statistics.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/odl/util/statistics.py b/odl/util/statistics.py index f90df7e7009..48264ac3f4c 100644 --- a/odl/util/statistics.py +++ b/odl/util/statistics.py @@ -30,16 +30,16 @@ def mse(true, noisy): """Return the Mean Squared Errror. - + Parameters ---------- true : array-like noisy : array-like - + Returns ------- - psnr : float - + mse : float + Examples -------- >>> true = [1, 1, 1, 1, 1] @@ -55,16 +55,16 @@ def mse(true, noisy): def psnr(true, noisy): """Return the Peak signal-to-noise ratio. - + Parameters ---------- true : array-like noisy : array-like - + Returns ------- psnr : float - + Examples -------- >>> true = [1, 1, 1, 1, 1] @@ -72,23 +72,23 @@ def psnr(true, noisy): >>> result = psnr(true, noisy) >>> print('{:.3f}'.format(result)) 6.990 - + If true == noisy, the result is infinite - + >>> psnr([1, 1], [1, 1]) inf - + If `true == 0` but `noisy != 0`, the result is negative infinity - + >>> psnr(0, 1) -inf """ true = np.asarray(true) noisy = np.asarray(noisy) - + mse_result = mse(true, noisy) max_true = np.max(np.abs(true)) - + if mse_result == 0: return np.inf elif max_true == 0: From c9ef813ef37b1008801508c3e16f920922cf9100 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 13 Apr 2017 18:53:37 +0200 Subject: [PATCH 50/70] MAINT: Cleanup tensorflow examples --- .../tensorflow/tensorflow_layer_matrix.py | 42 +++--- .../tensorflow_layer_ray_transform.py | 48 ++++--- .../tensorflow/tensorflow_tomography.py | 54 ++++---- .../tensorflow/tensorflow_tomography_cnn.py | 61 -------- .../tensorflow_tomography_cnn_gradient.py | 130 ------------------ 5 files changed, 84 insertions(+), 251 deletions(-) delete mode 100644 examples/contrib/tensorflow/tensorflow_tomography_cnn.py delete mode 100644 examples/contrib/tensorflow/tensorflow_tomography_cnn_gradient.py diff --git a/examples/contrib/tensorflow/tensorflow_layer_matrix.py b/examples/contrib/tensorflow/tensorflow_layer_matrix.py index 7780230a072..41f32998f13 100644 --- a/examples/contrib/tensorflow/tensorflow_layer_matrix.py +++ b/examples/contrib/tensorflow/tensorflow_layer_matrix.py @@ -4,29 +4,35 @@ import numpy as np import odl +sess = tf.InteractiveSession() -with tf.Session() as sess: - tf.global_variables_initializer().run() +# Define ODL operator +matrix = np.array([[1, 2], + [0, 0], + [0, 1]], dtype='float32') +odl_op = odl.MatrixOperator(matrix) - matrix = np.array([[1, 2], - [0, 0], - [0, 1]], dtype='float32') - odl_op = odl.MatrixOperator(matrix) +# Define evaluation points +x = [1., 2.] +z = [1., 2., 3.] - x = tf.constant([1., 2.]) - z = tf.constant([1., 2., 3.]) +# Add empty axes for batch and channel +x_tf = tf.constant(x)[None, ..., None] +z_tf = tf.constant(z)[None, ..., None] - odl_op_layer = odl.as_tensorflow_layer(odl_op, 'MatrixOperator') - y = odl_op_layer(x) +# Create tensorflow layer from odl operator +odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( + odl_op, 'MatrixOperator') +y_tf = odl_op_layer(x_tf) - # Evaluate using tensorflow - print(y.eval()) +# Evaluate using tensorflow +print(y_tf.eval().ravel()) - # Compare result with pure ODL - print(odl_op(x.eval())) +# Compare result with pure ODL +print(odl_op(x)) - # Evaluate the adjoint of the derivative, called gradient in tensorflow - print(tf.gradients(y, [x], z)[0].eval()) +# Evaluate the adjoint of the derivative, called gradient in tensorflow +print(tf.gradients(y_tf, [x_tf], z_tf)[0].eval().ravel()) - # Compare result with pure ODL - print(odl_op.derivative(x.eval()).adjoint(z.eval())) +# Compare result with pure ODL +print(odl_op.derivative(x).adjoint(z)) diff --git a/examples/contrib/tensorflow/tensorflow_layer_ray_transform.py b/examples/contrib/tensorflow/tensorflow_layer_ray_transform.py index 66ad0e2124c..3197e214200 100644 --- a/examples/contrib/tensorflow/tensorflow_layer_ray_transform.py +++ b/examples/contrib/tensorflow/tensorflow_layer_ray_transform.py @@ -4,29 +4,41 @@ import numpy as np import odl +sess = tf.InteractiveSession() -with tf.Session() as sess: - tf.global_variables_initializer().run() +tf.global_variables_initializer().run() - space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], - dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space) - ray_transform = odl.tomo.RayTransform(space, geometry) +space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], + dtype='float32') +geometry = odl.tomo.parallel_beam_geometry(space) +ray_transform = odl.tomo.RayTransform(space, geometry) - x = tf.constant(np.asarray(ray_transform.domain.one())) - z = tf.constant(np.asarray(ray_transform.range.one())) +x = tf.constant(np.asarray(ray_transform.domain.one())) +z = tf.constant(np.asarray(ray_transform.range.one())) - odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') - y = odl_op_layer(x) +# Create tensorflow layer from odl operator +odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( + ray_transform, 'RayTransform') - # Evaluate using tensorflow - print(y.eval()) - # Compare result with pure ODL - print(ray_transform(x.eval())) +# Add empty axes for batch and channel +x = x[None, ..., None] +z = z[None, ..., None] - # Evaluate the adjoint of the derivative, called gradient in tensorflow - print(tf.gradients(y, [x], z)[0].eval()) +# Lazily apply operator in tensorflow +y = odl_op_layer(x) - # Compare result with pure ODL - print(ray_transform.derivative(x.eval()).adjoint(z.eval())) +# Evaluate using tensorflow +print(y.eval()) + +# Compare result with pure ODL +print(ray_transform(x.eval())) + +# Evaluate the adjoint of the derivative, called gradient in tensorflow +# We need to scale by cell size to get correct value since the derivative +# in tensorflow uses unweighted spaces. +scale = ray_transform.range.cell_volume / ray_transform.domain.cell_volume +print(tf.gradients(y, [x], z)[0].eval() * scale) + +# Compare result with pure ODL +print(ray_transform.derivative(x.eval()).adjoint(z.eval())) diff --git a/examples/contrib/tensorflow/tensorflow_tomography.py b/examples/contrib/tensorflow/tensorflow_tomography.py index bea55677564..23a377c2d14 100644 --- a/examples/contrib/tensorflow/tensorflow_tomography.py +++ b/examples/contrib/tensorflow/tensorflow_tomography.py @@ -5,35 +5,41 @@ import odl -with tf.Session() as sess: - # Create ODL data structures - space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], - dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space) - ray_transform = odl.tomo.RayTransform(space, geometry) - phantom = odl.phantom.shepp_logan(space, True) - data = ray_transform(phantom) +sess = tf.InteractiveSession() - # Create tensorflow layer from odl operator - odl_op_layer = odl.as_tensorflow_layer(ray_transform, 'RayTransform') - x = tf.Variable(tf.constant(0.0, shape=space.shape), name="x") +# Create ODL data structures +space = odl.uniform_discr([-64, -64], [64, 64], [128, 128], + dtype='float32') +geometry = odl.tomo.parallel_beam_geometry(space) +ray_transform = odl.tomo.RayTransform(space, geometry) +phantom = odl.phantom.shepp_logan(space, True) +data = ray_transform(phantom) - # Create constant right hand side - y = tf.constant(np.asarray(data)) +# Create tensorflow layer from odl operator +odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( + ray_transform, 'RayTransform') +x = tf.Variable(tf.constant(0.0, shape=space.shape), name="x") - # Define l2 loss function - loss = tf.nn.l2_loss(odl_op_layer(x) - y) +# Create constant right hand side +y = tf.constant(np.asarray(data)) - # Train using the adam optimizer - optimizer = tf.train.AdamOptimizer(1e-1).minimize(loss) +# Add empty axes for batch and channel +x = x[None, ..., None] +y = y[None, ..., None] - # Initialize all TF variables - tf.global_variables_initializer().run() +# Define l2 loss function +loss = tf.nn.l2_loss(odl_op_layer(x) - y) - # Solve with an ODL callback to see what happens - callback = odl.solvers.CallbackShow() +# Train using the adam optimizer +optimizer = tf.train.AdamOptimizer(1e-1).minimize(loss) - for i in range(100): - sess.run(optimizer) +# Initialize all TF variables +tf.global_variables_initializer().run() - callback(space.element(x.eval())) +# Solve with an ODL callback to see what happens +callback = odl.solvers.CallbackShow() + +for i in range(100): + sess.run(optimizer) + + callback(space.element(x.eval())) diff --git a/examples/contrib/tensorflow/tensorflow_tomography_cnn.py b/examples/contrib/tensorflow/tensorflow_tomography_cnn.py deleted file mode 100644 index 2e52ef29d9d..00000000000 --- a/examples/contrib/tensorflow/tensorflow_tomography_cnn.py +++ /dev/null @@ -1,61 +0,0 @@ -"""Example of how to solve a simple tomography problem using tensorflow.""" - -import tensorflow as tf -import numpy as np -import odl - - -def conv2d(x, W): - return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME') - - -with tf.Session() as sess: - # Create ODL data structures - space = odl.uniform_discr([-64, -64], [64, 64], [512, 512], - dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space, angles=20) - ray_transform = odl.tomo.RayTransform(space, geometry) - phantom = odl.phantom.shepp_logan(space, True) - data = ray_transform(phantom) - fbp = odl.tomo.fbp_op(ray_transform)(data) - - # Create tensorflow layer from odl operator - odl_op_layer = odl.as_tensorflow_layer(ray_transform, - 'RayTransform') - x = tf.constant(np.asarray(fbp), name="x") - x = x[None, :, :, None] - - # Create constant right hand side - y = tf.constant(np.asarray(data)) - - init = np.random.randn(*[3, 3, 1, 32]) * 0.1 - W = tf.Variable(tf.constant(init, dtype='float32')) - x = tf.nn.relu(conv2d(x, W)) - - for i in range(2): - init = np.random.randn(*[3, 3, 32, 32]) * 0.1 - W = tf.Variable(tf.constant(init, dtype='float32')) - x = tf.nn.relu(conv2d(x, W)) - - init = np.random.randn(*[3, 3, 32, 1]) * 0.1 - W = tf.Variable(tf.constant(init, dtype='float32')) - x = tf.nn.relu(conv2d(x, W)) - - # Reshape for ODL - x = x[0, :, :, 0] - - # Define l2 loss function - loss = tf.nn.l2_loss(odl_op_layer(x) - y) - - # Train using the adam optimizer - optimizer = tf.train.AdamOptimizer(1e-3).minimize(loss) - - # Initialize all TF variables - tf.global_variables_initializer().run() - - # Solve with an ODL callback to see what happens - callback = odl.solvers.CallbackShow(clim=[0, 1]) - - for i in range(1000): - sess.run(optimizer) - callback(space.element(x.eval())) diff --git a/examples/contrib/tensorflow/tensorflow_tomography_cnn_gradient.py b/examples/contrib/tensorflow/tensorflow_tomography_cnn_gradient.py deleted file mode 100644 index dcc314f1e2b..00000000000 --- a/examples/contrib/tensorflow/tensorflow_tomography_cnn_gradient.py +++ /dev/null @@ -1,130 +0,0 @@ -"""Example of how to solve a simple tomography problem using tensorflow.""" - -import tensorflow as tf -import numpy as np -import odl - - -def random_ellipse(): - return (np.random.rand() - 0.3, - np.random.exponential() * 0.2, np.random.exponential() * 0.2, - np.random.rand() - 0.5, np.random.rand() - 0.5, - np.random.rand() * 2 * np.pi) - - -def random_phantom(spc): - n = np.random.poisson(10) - ellipses = [random_ellipse() for _ in range(n)] - return odl.phantom.ellipsoid_phantom(spc, ellipses) - - -def conv2d(x, W): - return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME') - - -def var(x): - return tf.Variable(tf.constant(x, dtype='float32')) - - -def create_variable(name, shape, stddev=0.01): - variable = tf.Variable(tf.truncated_normal(shape, stddev=stddev), name=name) - return variable - - -with tf.Session() as sess: - - # Create ODL data structures - size = 128 - space = odl.uniform_discr([-64, -64], [64, 64], [size, size], - dtype='float32') - geometry = odl.tomo.parallel_beam_geometry(space) - ray_transform = odl.tomo.RayTransform(space, geometry) - - - # Create tensorflow layer from odl operator - odl_op_layer = odl.as_tensorflow_layer(ray_transform, - 'RayTransform') - odl_op_layer_adjoint = odl.as_tensorflow_layer(ray_transform.adjoint, - 'RayTransformAdjoint') - - n_data = 50 - x_arr = np.empty((n_data, space.shape[0], space.shape[1], 1), dtype='float32') - y_arr = np.empty((n_data, ray_transform.range.shape[0], ray_transform.range.shape[1], 1), dtype='float32') - x_true_arr = np.empty((n_data, space.shape[0], space.shape[1], 1), dtype='float32') - - for i in range(n_data): - if i == 0: - phantom = odl.phantom.shepp_logan(space, True) - else: - phantom = random_phantom(space) - data = ray_transform(phantom) - noisy_data = data + odl.phantom.white_noise(ray_transform.range) * np.mean(data) * 0.1 - fbp = odl.tomo.fbp_op(ray_transform)(noisy_data) - - x_arr[i] = np.asarray(fbp)[..., None] - x_true_arr[i] = np.asarray(phantom)[..., None] - y_arr[i] = np.asarray(noisy_data)[..., None] - - x_0 = tf.placeholder(tf.float32, shape=[None, size, size, 1], name="x_0") - x_true = tf.placeholder(tf.float32, shape=[None, size, size, 1], name="x_true") - y = tf.placeholder(tf.float32, shape=[None, ray_transform.range.shape[0], ray_transform.range.shape[1], 1], name="y") - - s = tf.fill([tf.shape(x_0)[0], size, size, 5], np.float32(0.0), name="s") - - # Create constant right hand side - - w1 = tf.Variable(tf.truncated_normal([3, 3, 7, 32], stddev=0.01)) - b1 = var(np.ones(32) * 0.1) - - w2 = tf.Variable(tf.truncated_normal([3, 3, 32, 32], stddev=0.01)) - b2 = var(np.ones(32) * 0.1) - - w3 = tf.Variable(tf.truncated_normal([3, 3, 32, 6], stddev=0.01)) - b3 = var(np.random.randn(6) * 0.001) - - n_iter = 5 - - x = x_0 - loss = 0 - for i in range(n_iter): - gradx = odl_op_layer_adjoint(odl_op_layer(x) - y) - - update = tf.concat([x, gradx, s], axis=3) - - update = tf.nn.relu(conv2d(update, w1) + b1) - update = tf.nn.relu(conv2d(update, w2) + b2) - update = tf.nn.dropout(update, 0.8) - update = tf.nn.tanh(conv2d(update, w3) + b3) - - ds = update[..., 1:] - dx = update[..., 0][..., None] - - s = s + ds - x = x + dx - - loss = loss + tf.nn.l2_loss(x - x_true) * (2. ** (i - n_iter)) / n_data - - # Train using the adam optimizer - learning_rate = tf.placeholder(tf.float32, shape=[]) - optimizer = tf.train.AdamOptimizer(learning_rate).minimize(loss) - - # Initialize all TF variables - tf.global_variables_initializer().run() - - # Solve with an ODL callback to see what happens - callback = odl.solvers.CallbackShow() - - for i in range(10000): - _, loss_training = sess.run([optimizer, loss], - feed_dict={learning_rate: 0.001, - x_0: x_arr[1:], - x_true: x_true_arr[1:], - y: y_arr[1:]}) - - x_value, loss_value = sess.run([x, loss], - feed_dict={x_0: x_arr[0:1], - x_true: x_true_arr[0:1], - y: y_arr[0:1]}) - - print('iter={}, training loss={}, validation loss={}'.format(i, loss_training, loss_value)) - callback(space.element(x_value[0])) From 1fe5ab244e2b86b26666d4a5e4dedef793bfbc7b Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Wed, 19 Apr 2017 23:38:38 +0200 Subject: [PATCH 51/70] API: Change shape of productspace and add asarray --- .../tensorflow_layer_productspace.py | 38 +++++++++++++++++++ odl/space/pspace.py | 33 +++++++++++++++- 2 files changed, 69 insertions(+), 2 deletions(-) create mode 100644 examples/contrib/tensorflow/tensorflow_layer_productspace.py diff --git a/examples/contrib/tensorflow/tensorflow_layer_productspace.py b/examples/contrib/tensorflow/tensorflow_layer_productspace.py new file mode 100644 index 00000000000..7664ed619bd --- /dev/null +++ b/examples/contrib/tensorflow/tensorflow_layer_productspace.py @@ -0,0 +1,38 @@ +"""Example of how to convert an ODL operator to a tensorflow layer.""" + +import tensorflow as tf +import numpy as np +import odl + +sess = tf.InteractiveSession() + +# Define ODL operator +space = odl.uniform_discr([0, 0], [1, 1], [10, 10], dtype='float32') + +odl_op = odl.Gradient(space) + +# Define evaluation points +x = odl_op.domain.one() +z = odl_op.range.one() + +# Add empty axes for batch and channel +x_tf = tf.ones([1, 10, 10, 1]) +z_tf = tf.ones([1, 2, 10, 10, 1]) + +# Create tensorflow layer from odl operator +odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( + odl_op, 'Gradient') +y_tf = odl_op_layer(x_tf) + +# Evaluate using tensorflow +print(y_tf.eval().squeeze()) + +# Compare result with pure ODL +print(odl_op(x)) + +# Evaluate the adjoint of the derivative, called gradient in tensorflow +scale = 1 / space.cell_volume +print(tf.gradients(y_tf, [x_tf], z_tf)[0].eval().squeeze() * scale) + +# Compare result with pure ODL +print(odl_op.derivative(x).adjoint(z)) diff --git a/odl/space/pspace.py b/odl/space/pspace.py index 9ee2f85a6f5..4dc5d7a7aed 100644 --- a/odl/space/pspace.py +++ b/odl/space/pspace.py @@ -309,7 +309,10 @@ def __len__(self): def shape(self): """Number of spaces per axis.""" # Currently supporting only 1d product spaces - return (self.size,) + if not self.is_power_space: + return NotImplementedError('shape only defined for power spaces') + else: + return (self.size,) + self[0].shape @property def spaces(self): @@ -547,7 +550,7 @@ def __eq__(self, other): return True else: return (isinstance(other, ProductSpace) and - self.shape == other.shape and + len(self) == len(other) and self.weighting == other.weighting and all(x == y for x, y in zip(self.spaces, other.spaces))) @@ -626,6 +629,11 @@ def parts(self): """Parts of this product space element.""" return self.__parts + @property + def shape(self): + """Number of spaces per axis.""" + return self.space.shape + @property def size(self): """Number of factors of this element's space.""" @@ -705,6 +713,27 @@ def __setitem__(self, indices, values): for p, v in zip(indexed_parts, values): p[:] = v + def __array__(self): + """An array representation of ``self``. + + Only available if `is_power_space` is True. + + Examples + -------- + >>> spc = odl.ProductSpace(odl.rn(2), 2) + >>> x = spc.one() + >>> np.asarray(x) + array([[ 1., 1.], + [ 1., 1.]]) + """ + if not self.space.is_power_space: + return NotImplemented + else: + arr = np.zeros(self.shape, self.dtype) + for i in range(len(self)): + arr[i] = np.asarray(self[i]) + return arr + @property def ufuncs(self): """`ProductSpaceUfuncs`, access to Numpy style ufuncs. From ba647fde239aad1e8a9f2ad22a608d4e8382630c Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 20 Apr 2017 14:38:48 +0200 Subject: [PATCH 52/70] ENH: Add estimate_noise --- odl/util/statistics.py | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/odl/util/statistics.py b/odl/util/statistics.py index 48264ac3f4c..db04c55506c 100644 --- a/odl/util/statistics.py +++ b/odl/util/statistics.py @@ -25,7 +25,7 @@ import numpy as np -__all__ = ('psnr',) +__all__ = ('mse', 'psnr', 'estimate_noise') def mse(true, noisy): @@ -97,6 +97,43 @@ def psnr(true, noisy): return 20 * np.log10(max_true) - 10 * np.log10(mse_result) +def estimate_noise(img): + """Estimate noise level in image. + + Parameters + ---------- + img : array-lik + + Returns + ------- + noise : float + estimate of standard deviation in image + + Examples + -------- + Create image with noise 1.0, verify result + + >>> img = np.random.randn(10, 10) + >>> noise = estimate_noise(img) # should be about 1 + + See Also + -------- + Reference: J. Immerkaer, “Fast Noise Variance Estimation” + """ + import scipy.signal + import functools + img = np.asarray(img, dtype='float') + + M = functools.reduce(np.add.outer, [[-1, 2, -1]] * img.ndim) + + convolved = scipy.signal.convolve(img, M, 'valid') + conv_var = np.sum(np.square(convolved)) + + sigma = np.sqrt(conv_var / (np.sum(M**2) * np.prod(convolved.shape))) + + return sigma + + if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From 188d9dd97bb08eda1582ad8ab859b7ce8ad2c600 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 12 May 2017 16:18:34 +0200 Subject: [PATCH 53/70] ENH: Add option to normalize PSNR --- odl/util/statistics.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/odl/util/statistics.py b/odl/util/statistics.py index db04c55506c..f9d15726021 100644 --- a/odl/util/statistics.py +++ b/odl/util/statistics.py @@ -53,13 +53,16 @@ def mse(true, noisy): return np.mean((true - noisy) ** 2) -def psnr(true, noisy): +def psnr(true, noisy, normalize=False): """Return the Peak signal-to-noise ratio. Parameters ---------- true : array-like noisy : array-like + normalize : bool + If true, normalizes the true and noisy signal to have the same mean + and variance before comparison. Returns ------- @@ -86,6 +89,10 @@ def psnr(true, noisy): true = np.asarray(true) noisy = np.asarray(noisy) + if normalize: + noisy = noisy + np.mean(true - noisy) + noisy = noisy * np.std(true) / np.std(noisy) + mse_result = mse(true, noisy) max_true = np.max(np.abs(true)) From 34f11fc3f7cd8b838fa0a6abe72f72b324f9dab3 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Mon, 22 May 2017 18:12:32 +0200 Subject: [PATCH 54/70] ENH: Add ADAM solver, see #984 --- odl/solvers/smooth/gradient.py | 73 +++++++++++++++++++- odl/test/solvers/iterative/iterative_test.py | 7 ++ odl/test/solvers/smooth/smooth_test.py | 9 +++ 3 files changed, 88 insertions(+), 1 deletion(-) diff --git a/odl/solvers/smooth/gradient.py b/odl/solvers/smooth/gradient.py index 1b558b71ddd..aebbf065f67 100644 --- a/odl/solvers/smooth/gradient.py +++ b/odl/solvers/smooth/gradient.py @@ -18,7 +18,7 @@ from odl.solvers.util import ConstantLineSearch -__all__ = ('steepest_descent',) +__all__ = ('steepest_descent', 'adam') # TODO: update all docs @@ -109,6 +109,77 @@ def steepest_descent(f, x, line_search=1.0, maxiter=1000, tol=1e-16, callback(x) +def adam(f, x, learning_rate=1e-3, beta1=0.9, beta2=0.999, eps=1e-8, + maxiter=1000, tol=1e-16, callback=None): + """ADAM method to minimize an objective function. + + General implementation of ADAM for solving + + .. math:: + \min f(x) + + The algorithm is intended for unconstrained problems. + + The algorithm is described in + `Adam: A Method for Stochastic Optimization + `_. All parameter names are taken from + that article. + + Parameters + ---------- + f : `Functional` + Goal functional. Needs to have ``f.gradient``. + x : ``f.domain`` element + Starting point of the iteration + learning_rate : float, optional + Step length of the method. + beta1 : float, optional + Update rate for first order moment estimate. + beta2 : float, optional + Update rate for second order moment estimate. + eps : float, optional + A small constant for numerical stability. + maxiter : int, optional + Maximum number of iterations. + tol : float, optional + Tolerance that should be used for terminating the iteration. + callback : callable, optional + Object executing code per iteration, e.g. plotting each iterate + + See Also + -------- + odl.solvers.smooth.gradient.steepest_descent : simple steepest descent + odl.solvers.iterative.iterative.landweber : + Optimized solver for the case ``f(x) = ||Ax - b||_2^2`` + odl.solvers.iterative.iterative.conjugate_gradient : + Optimized solver for the case ``f(x) = x^T Ax - 2 x^T b`` + """ + grad = f.gradient + if x not in grad.domain: + raise TypeError('`x` {!r} is not in the domain of `grad` {!r}' + ''.format(x, grad.domain)) + + m = grad.domain.zero() + v = grad.domain.zero() + + grad_x = grad.range.element() + for _ in range(maxiter): + grad(x, out=grad_x) + + if grad_x.norm() < tol: + return + + m.lincomb(beta1, m, 1 - beta1, grad_x) + v.lincomb(beta2, v, 1 - beta2, grad_x ** 2) + + step = learning_rate * np.sqrt(1 - beta2) / (1 - beta1) + + x.lincomb(1, x, -step, m / np.sqrt(v + eps)) + + if callback is not None: + callback(x) + + if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests diff --git a/odl/test/solvers/iterative/iterative_test.py b/odl/test/solvers/iterative/iterative_test.py index 1b1aed21e7b..edab585943f 100644 --- a/odl/test/solvers/iterative/iterative_test.py +++ b/odl/test/solvers/iterative/iterative_test.py @@ -18,6 +18,7 @@ # Find the valid projectors @pytest.fixture(scope="module", params=['steepest_descent', + 'adam', 'landweber', 'conjugate_gradient', 'conjugate_gradient_normal', @@ -34,6 +35,12 @@ def solver(op, x, rhs): func = odl.solvers.L2NormSquared(op.domain) * (op - rhs) odl.solvers.steepest_descent(func, x, line_search=0.5 / norm2) + elif solver_name == 'adam': + def solver(op, x, rhs): + norm2 = op.adjoint(op(x)).norm() / x.norm() + func = odl.solvers.L2NormSquared(op.domain) * (op - rhs) + + odl.solvers.adam(func, x, learning_rate=4.0 / norm2, maxiter=150) elif solver_name == 'landweber': def solver(op, x, rhs): norm2 = op.adjoint(op(x)).norm() / x.norm() diff --git a/odl/test/solvers/smooth/smooth_test.py b/odl/test/solvers/smooth/smooth_test.py index 40a4cfa49b8..9c21f055e2d 100644 --- a/odl/test/solvers/smooth/smooth_test.py +++ b/odl/test/solvers/smooth/smooth_test.py @@ -123,6 +123,15 @@ def test_steepest_descent(functional): assert functional(x) < 1e-3 +def test_adam(functional): + """Test the ``adam`` solver.""" + + x = functional.domain.one() + odl.solvers.adam(functional, x, tol=1e-2, learning_rate=0.5) + + assert functional(x) < 1e-3 + + def test_conjguate_gradient_nonlinear(functional, nonlinear_cg_beta): """Test the ``conjugate_gradient_nonlinear`` solver.""" line_search = odl.solvers.BacktrackingLineSearch(functional) From 7722357968b3c222a9261bb760b2a9979d34793c Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 1 Aug 2017 13:29:12 +0200 Subject: [PATCH 55/70] MAINT: Clean up tensorflow package structure --- odl/contrib/tensorflow/__init__.py | 20 +++ .../examples}/tensorflow_layer_matrix.py | 0 .../tensorflow_layer_productspace.py | 0 .../tensorflow_layer_ray_transform.py | 0 .../examples/tensorflow_operator.py | 142 +++++++++++++++++ .../examples}/tensorflow_tomography.py | 0 .../{tensorflow.py => tensorflow/layer.py} | 134 +--------------- odl/contrib/tensorflow/space.py | 149 ++++++++++++++++++ 8 files changed, 315 insertions(+), 130 deletions(-) create mode 100644 odl/contrib/tensorflow/__init__.py rename {examples/contrib/tensorflow => odl/contrib/tensorflow/examples}/tensorflow_layer_matrix.py (100%) rename {examples/contrib/tensorflow => odl/contrib/tensorflow/examples}/tensorflow_layer_productspace.py (100%) rename {examples/contrib/tensorflow => odl/contrib/tensorflow/examples}/tensorflow_layer_ray_transform.py (100%) create mode 100644 odl/contrib/tensorflow/examples/tensorflow_operator.py rename {examples/contrib/tensorflow => odl/contrib/tensorflow/examples}/tensorflow_tomography.py (100%) rename odl/contrib/{tensorflow.py => tensorflow/layer.py} (76%) create mode 100644 odl/contrib/tensorflow/space.py diff --git a/odl/contrib/tensorflow/__init__.py b/odl/contrib/tensorflow/__init__.py new file mode 100644 index 00000000000..3cfe4655e58 --- /dev/null +++ b/odl/contrib/tensorflow/__init__.py @@ -0,0 +1,20 @@ +# Copyright 2014-2017 The ODL contributors +# +# This file is part of ODL. +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. + +"""ODL interactions with tensorflow.""" + + +from __future__ import absolute_import + +__all__ = () + +from .layer import * +__all__ += layer. + +from .space import * +__all__ += space.__all__ diff --git a/examples/contrib/tensorflow/tensorflow_layer_matrix.py b/odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py similarity index 100% rename from examples/contrib/tensorflow/tensorflow_layer_matrix.py rename to odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py diff --git a/examples/contrib/tensorflow/tensorflow_layer_productspace.py b/odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py similarity index 100% rename from examples/contrib/tensorflow/tensorflow_layer_productspace.py rename to odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py diff --git a/examples/contrib/tensorflow/tensorflow_layer_ray_transform.py b/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py similarity index 100% rename from examples/contrib/tensorflow/tensorflow_layer_ray_transform.py rename to odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py diff --git a/odl/contrib/tensorflow/examples/tensorflow_operator.py b/odl/contrib/tensorflow/examples/tensorflow_operator.py new file mode 100644 index 00000000000..5ca9e8aa855 --- /dev/null +++ b/odl/contrib/tensorflow/examples/tensorflow_operator.py @@ -0,0 +1,142 @@ +"""Example of wrapping tensorflow to create an ODL operator.""" +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.legend_handler import HandlerLine2D +import tensorflow as tf +import odl + + +sess = tf.InteractiveSession() + +def G(x): + if(np.array(x.shape).size==1): + return 1./(np.pi**.5)*np.exp(-x**2/2.) + elif(np.array(x.shape).size==2): + return 1./(np.pi**.5)*np.exp(-np.sum(x**2,axis=1)/2.) + else: + print('Error: x should be an MxN array') + return + +def G_sigma_m_x_m(x,sigma_m,x_m): + return G((x-x_m)/sigma_m) + +def fp(Arg,x): + # Arg = 4xM array + # Arg = [[alpha_j],[sigma_j],[x1_j],[x2_j]] + fp_temp = np.zeros(x.shape[0]) + M = Arg.shape[1] + alpha_j = Arg[0] + sigma_j = Arg[1] + x_j = Arg[2:].T + + for m in np.arange(M): + fp_temp = fp_temp + alpha_j[m]*G((x-x_j[m])/sigma_j[m]) + + return fp_temp + + +def R_G(p): + return tf.exp(-p ** 2. / 2.) + + +def R_G_sigma_m_x_m(theta, r, sigma_m, x_m, y_m): + return R_G((r - (tf.cos(theta) * x_m + tf.sin(theta) * y_m)) / sigma_m) + + +def R_fp(alpha, sigma, x, y, theta, r): + R_fp_temp = R_G_sigma_m_x_m(theta[:, None], r[:, None], + sigma[None, :], x[None, :], y[None, :]) + + R_fp_temp = alpha[None, :] * R_fp_temp + + return tf.reduce_sum(R_fp_temp, axis=-1) + + +class GaussianRayTransform(odl.Operator): + def __init__(self, domain, range, theta, p): + size = domain[0].size + self.alpha_ph = tf.placeholder(tf.float32, shape=[size]) + self.sigma_ph = tf.placeholder(tf.float32, shape=[size]) + self.x_ph = tf.placeholder(tf.float32, shape=[size]) + self.y_ph = tf.placeholder(tf.float32, shape=[size]) + theta = tf.constant(theta, dtype=tf.float32) + p = tf.constant(p, dtype=tf.float32) + + self.y = tf.placeholder(tf.float32, shape=range.shape) + + self.result = R_fp(self.alpha_ph, self.sigma_ph, self.x_ph, self.y_ph, theta, p) + self.gradient = tf.gradients(self.result, + [self.alpha_ph, self.sigma_ph, self.x_ph, self.y_ph], + [self.y]) + self.gradient = [gi * range.cell_volume for gi in self.gradient] + + odl.Operator.__init__(self, domain, range) + + def _call(self, x): + result = sess.run(self.result, + feed_dict={self.alpha_ph: np.asarray(x[0]), + self.sigma_ph: np.asarray(x[1]), + self.x_ph: np.asarray(x[2]), + self.y_ph: np.asarray(x[3])}) + + return result.reshape(self.range.shape).T + + def derivative(self, x): + op = self + + class GaussianRayTransformDerivative(odl.Operator): + @property + def adjoint(self): + class GaussianRayTransformDerivativeAdjoint(odl.Operator): + def _call(self, y): + result = sess.run(op.gradient, + feed_dict={op.alpha_ph: np.asarray(x[0]), + op.sigma_ph: np.asarray(x[1]), + op.x_ph: np.asarray(x[2]), + op.y_ph: np.asarray(x[3]), + op.y: np.asarray(y).T}) + + return result + + return GaussianRayTransformDerivativeAdjoint(self.range, + self.domain, + linear=True) + + return GaussianRayTransformDerivative(self.domain, self.range, + linear=True) + + +if __name__ == '__main__': + n = 201 + x,y = np.meshgrid(np.linspace(-2,2,n)*5,np.linspace(-2,2,n)*5) + X = np.row_stack((x.flatten(),y.flatten())).T + + M = 13 + r = np.linspace(0,2,M)*4 + x_theta = np.linspace(0,2.*np.pi,M) + x_ini = r*np.array([np.cos(x_theta),np.sin(x_theta)]) + #plt.plot(x_ini[0],x_ini[1],'.') + #plt.show() + alpha_ini = np.ones(M) + sigma_ini = np.ones(M) + Arg_ini = np.row_stack((alpha_ini,sigma_ini,x_ini)) + + theta,p = np.meshgrid(np.linspace(0,2.*np.pi,n),np.linspace(-20,20,n)) + Theta = theta.flatten() + P = p.flatten() + + import odl + rn = odl.rn(M) + domain = rn ** 4 + + ran = odl.uniform_discr([0, -20], [2 * np.pi, 20], [n, n]) + operator = GaussianRayTransform(domain, ran, Theta, P) + + x = domain.element([alpha_ini, sigma_ini, x_ini[0], x_ini[1]]) + + y = operator(x) + y.show() + + deriv = operator.derivative(x) + + adjy = deriv.adjoint(y) diff --git a/examples/contrib/tensorflow/tensorflow_tomography.py b/odl/contrib/tensorflow/examples/tensorflow_tomography.py similarity index 100% rename from examples/contrib/tensorflow/tensorflow_tomography.py rename to odl/contrib/tensorflow/examples/tensorflow_tomography.py diff --git a/odl/contrib/tensorflow.py b/odl/contrib/tensorflow/layer.py similarity index 76% rename from odl/contrib/tensorflow.py rename to odl/contrib/tensorflow/layer.py index addac92bed8..64e7640c2a0 100644 --- a/odl/contrib/tensorflow.py +++ b/odl/contrib/tensorflow/layer.py @@ -15,7 +15,7 @@ # You should have received a copy of the GNU General Public License # along with ODL. If not, see . -"""Utilities for interfacing ODL with Tensorflow.""" +"""Utilities for converting ODL operators to tensorflow layers.""" # Imports for common Python 2/3 codebase from __future__ import print_function, division, absolute_import @@ -23,20 +23,13 @@ from future import standard_library standard_library.install_aliases() -from odl.set import LinearSpace, LinearSpaceElement, RealNumbers -from odl.operator import Operator - import numpy as np import uuid -try: - # TODO: Lazy import to improve performance - import tensorflow as tf - TENSORFLOW_AVAILABLE = True -except: - TENSORFLOW_AVAILABLE = False +import tensorflow as tf +from tensorflow.python.framework import ops -__all__ = ('as_tensorflow_layer', 'TensorflowSpace', 'TensorflowSpaceOperator') +__all__ = ('as_tensorflow_layer') def as_tensorflow_layer(odl_op, default_name='ODLOperator', @@ -70,13 +63,7 @@ def as_tensorflow_layer(odl_op, default_name='ODLOperator', Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. The `dtype` of the tensor is the same as the respective ODL spaces. - """ - if not TENSORFLOW_AVAILABLE: - raise ImportError('Tensorflow not installed') - - from tensorflow.python.framework import ops - def py_func(func, inp, Tout, stateful=True, name=None, grad=None): """Define custom py_func which takes also a grad op as argument. @@ -393,119 +380,6 @@ def tensorflow_layer_grad(op, grad): return tensorflow_layer -class TensorflowSpace(LinearSpace): - - """A space of tensorflow Tensors.""" - - def __init__(self, shape, name='ODLTensorflowSpace'): - LinearSpace.__init__(self, RealNumbers()) - self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) - else si - for si in shape) - self.init_shape = tuple(si if si.value is not None - else tf.Dimension(1) - for si in self.shape) - self.name = name - - def _lincomb(self, a, x1, b, x2, out): - with tf.name_scope('{}_lincomb'.format(self.name)): - if x1 is x2: - # x1 is aligned with x2 -> out = (a+b)*x1 - out.data = (a + b) * x1.data - elif out is x1 and out is x2: - # All the vectors are aligned -> out = (a+b)*out - if (a + b) != 1: - out.data *= (a + b) - elif out is x1: - # out is aligned with x1 -> out = a*out + b*x2 - out.data = a * out.data + b * x2.data - elif out is x2: - # out is aligned with x2 -> out = a*x1 + b*out - out.data = a * x1.data + b * out.data - else: - # We have exhausted all alignment options, so x1 != x2 != out - # We now optimize for various values of a and b - if b == 0: - if a == 0: # Zero assignment -> out = 0 - out.data *= 0 - else: # Scaled copy -> out = a*x1 - out.data = a * x1.data - else: - if a == 0: # Scaled copy -> out = b*x2 - out.data = b * x2.data - elif a == 1: # No scaling in x1 -> out = x1 + b*x2 - out.data = x1.data + b * x2.data - else: # Generic case -> out = a*x1 + b*x2 - out.data = a * x1.data + b * x2.data - - def element(self, inp=None): - if inp in self: - return inp - elif inp is None: - return self.zero() - else: - return TensorflowSpaceElement(self, inp) - - def zero(self): - with tf.name_scope('{}_zero'.format(self.name)): - return self.element(tf.zeros(self.init_shape, - dtype=tf.float32)) - - def one(self): - with tf.name_scope('{}_one'.format(self.name)): - return self.element(tf.ones(self.init_shape, - dtype=tf.float32)) - - def __eq__(self, other): - return isinstance(other, TensorflowSpace) and other.shape == self.shape - - def __repr__(self): - return 'TensorflowSpace({})'.format(self.shape) - - -class TensorflowSpaceElement(LinearSpaceElement): - - """Elements in TensorflowSpace.""" - - def __init__(self, space, data): - LinearSpaceElement.__init__(self, space) - self.data = data - - @property - def data(self): - return self._data - - @data.setter - def data(self, value): - if isinstance(value, TensorflowSpaceElement): - raise Exception(value.data) - self._data = value - - def __repr__(self): - return '{}.element({})'.format(self.space, self.data) - - -class TensorflowSpaceOperator(Operator): - - """Wrap ODL operator so that it acts on TensorflowSpace elements.""" - - def __init__(self, domain, range, func, adjoint=None, linear=False): - Operator.__init__(self, domain, range, linear) - self.func = func - self.adjoint_func = adjoint - - def _call(self, x): - return self.func(x.data) - - @property - def adjoint(self): - return TensorflowSpaceOperator(self.range, - self.domain, - self.adjoint_func, - self.func, - self.is_linear) - - if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests diff --git a/odl/contrib/tensorflow/space.py b/odl/contrib/tensorflow/space.py new file mode 100644 index 00000000000..25bb9daf217 --- /dev/null +++ b/odl/contrib/tensorflow/space.py @@ -0,0 +1,149 @@ +# Copyright 2014-2017 The ODL development group +# +# This file is part of ODL. +# +# ODL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ODL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with ODL. If not, see . + +"""Utilities for converting ODL spaces to tensorflow layers.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import tensorflow as tf +from odl.set import LinearSpace, LinearSpaceElement, RealNumbers +from odl.operator import Operator + + +__all__ = ('TensorflowSpace', 'TensorflowSpaceOperator') + + +class TensorflowSpace(LinearSpace): + + """A space of tensorflow Tensors.""" + + def __init__(self, shape, name='ODLTensorflowSpace'): + LinearSpace.__init__(self, RealNumbers()) + self.shape = tuple(tf.Dimension(si) if not isinstance(si, tf.Dimension) + else si + for si in shape) + self.init_shape = tuple(si if si.value is not None + else tf.Dimension(1) + for si in self.shape) + self.name = name + + def _lincomb(self, a, x1, b, x2, out): + with tf.name_scope('{}_lincomb'.format(self.name)): + if x1 is x2: + # x1 is aligned with x2 -> out = (a+b)*x1 + out.data = (a + b) * x1.data + elif out is x1 and out is x2: + # All the vectors are aligned -> out = (a+b)*out + if (a + b) != 1: + out.data *= (a + b) + elif out is x1: + # out is aligned with x1 -> out = a*out + b*x2 + out.data = a * out.data + b * x2.data + elif out is x2: + # out is aligned with x2 -> out = a*x1 + b*out + out.data = a * x1.data + b * out.data + else: + # We have exhausted all alignment options, so x1 != x2 != out + # We now optimize for various values of a and b + if b == 0: + if a == 0: # Zero assignment -> out = 0 + out.data *= 0 + else: # Scaled copy -> out = a*x1 + out.data = a * x1.data + else: + if a == 0: # Scaled copy -> out = b*x2 + out.data = b * x2.data + elif a == 1: # No scaling in x1 -> out = x1 + b*x2 + out.data = x1.data + b * x2.data + else: # Generic case -> out = a*x1 + b*x2 + out.data = a * x1.data + b * x2.data + + def element(self, inp=None): + if inp in self: + return inp + elif inp is None: + return self.zero() + else: + return TensorflowSpaceElement(self, inp) + + def zero(self): + with tf.name_scope('{}_zero'.format(self.name)): + return self.element(tf.zeros(self.init_shape, + dtype=tf.float32)) + + def one(self): + with tf.name_scope('{}_one'.format(self.name)): + return self.element(tf.ones(self.init_shape, + dtype=tf.float32)) + + def __eq__(self, other): + return isinstance(other, TensorflowSpace) and other.shape == self.shape + + def __repr__(self): + return 'TensorflowSpace({})'.format(self.shape) + + +class TensorflowSpaceElement(LinearSpaceElement): + + """Elements in TensorflowSpace.""" + + def __init__(self, space, data): + LinearSpaceElement.__init__(self, space) + self.data = data + + @property + def data(self): + return self._data + + @data.setter + def data(self, value): + if isinstance(value, TensorflowSpaceElement): + raise Exception(value.data) + self._data = value + + def __repr__(self): + return '{}.element({})'.format(self.space, self.data) + + +class TensorflowSpaceOperator(Operator): + + """Wrap ODL operator so that it acts on TensorflowSpace elements.""" + + def __init__(self, domain, range, func, adjoint=None, linear=False): + Operator.__init__(self, domain, range, linear) + self.func = func + self.adjoint_func = adjoint + + def _call(self, x): + return self.func(x.data) + + @property + def adjoint(self): + return TensorflowSpaceOperator(self.range, + self.domain, + self.adjoint_func, + self.func, + self.is_linear) + + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() From f4028b88e6281304492334c22752cf0e7b53a361 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 1 Aug 2017 13:48:45 +0200 Subject: [PATCH 56/70] MAINT: fix bugs with tensorflow --- odl/contrib/tensorflow/__init__.py | 2 +- odl/contrib/tensorflow/layer.py | 2 +- odl/operator/oputils.py | 166 ----------------------------- 3 files changed, 2 insertions(+), 168 deletions(-) diff --git a/odl/contrib/tensorflow/__init__.py b/odl/contrib/tensorflow/__init__.py index 3cfe4655e58..5ae66c282a0 100644 --- a/odl/contrib/tensorflow/__init__.py +++ b/odl/contrib/tensorflow/__init__.py @@ -14,7 +14,7 @@ __all__ = () from .layer import * -__all__ += layer. +__all__ += layer.__all__ from .space import * __all__ += space.__all__ diff --git a/odl/contrib/tensorflow/layer.py b/odl/contrib/tensorflow/layer.py index 64e7640c2a0..4b49a576071 100644 --- a/odl/contrib/tensorflow/layer.py +++ b/odl/contrib/tensorflow/layer.py @@ -29,7 +29,7 @@ from tensorflow.python.framework import ops -__all__ = ('as_tensorflow_layer') +__all__ = ('as_tensorflow_layer',) def as_tensorflow_layer(odl_op, default_name='ODLOperator', diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 1582c3ce8b7..a16cde45e09 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -440,172 +440,6 @@ def adjoint(inp, out): norm_bound=norm_bound) -<<<<<<< b2414531c589a988c473ffc5cb7c9bebcb2b631c -def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): - """Convert ``Operator`` to tensorflow layer. - - Parameters - ---------- - odl_op : `Operator` - The operator that should be wrapped to a tensorflow layer. - name : str - Tensorflow name of the operator - differentiable : boolean - True if the operator should be differentiable, otherwise assumes that - the derivative is everywhere zero. - - Returns - ------- - tensorflow_layer : callable - Callable that, when called with an `tensorflow.Tensor` of shape - `(n, *odl_op.domain.shape, 1)` returns a tensor of shape - `(n, *odl_op.range.shape, 1)` where ``n`` is the number of batches. - Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. - The `dtype` of the tensor is the same as the respective ODL spaces. - - If ``odl_op`` implements `Operator.derivative` which in turn implements - `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and - gradients propagate as expected. - """ - import tensorflow as tf - from tensorflow.python.framework import ops - - def py_func(func, inp, Tout, stateful=True, name=None, grad=None): - """Define custom py_func which takes also a grad op as argument.""" - if grad is None: - return tf.py_func(func, inp, Tout, stateful=stateful, name=name) - else: - # Need to generate a unique name to avoid duplicates: - rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8)) - - tf.RegisterGradient(rnd_name)(grad) - g = tf.get_default_graph() - - if stateful: - override_name = 'PyFunc' - else: - override_name = 'PyFuncStateless' - - with g.gradient_override_map({override_name: rnd_name}): - return tf.py_func(func, inp, Tout, stateful=stateful, - name=name) - - def tensorflow_layer_grad_impl(x, dx): - """Implementation of the tensorflow gradient. - - Gradient in tensorflow is equivalent to the adjoint of the derivative - in ODL. - """ - x_shape = x.get_shape() - dx_shape = dx.get_shape() - try: - n_x = int(x_shape[0]) - fixed_size = True - except TypeError: - n_x = x_shape[0] - fixed_size = False - - in_shape = (n_x,) + odl_op.range.shape + (1,) - out_shape = (n_x,) + odl_op.domain.shape + (1,) - - # Validate input shape - assert x_shape[1:] == odl_op.domain.shape + (1,) - assert dx_shape[1:] == odl_op.range.shape + (1,) - - def _impl(x, dx): - if fixed_size: - x_out_shape = out_shape - assert x.shape == in_shape - assert dx.shape == out_shape - else: - x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == out_shape[1:] - assert dx.shape[1:] == in_shape[1:] - - out = np.empty(x_out_shape, odl_op.domain.dtype) - for i in range(x_out_shape[0]): - xi = x[i, ..., 0] - dxi = dx[i, ..., 0] - result = odl_op.derivative(xi).adjoint(dxi) - out[i] = np.asarray(result)[..., None] - - # TODO: Improve - scale = odl_op.domain.cell_volume / odl_op.range.cell_volume - out *= scale - return out - - with ops.name_scope(name + 'Grad', "ODLTensorflowLayerGrad", [x, dx]): - result = py_func(_impl, - [x, dx], - [tf.float32], - stateful=False) - - # We must manually set the output shape since tensorflow cannot - # figure it out - result = result[0] - result.set_shape(out_shape) - return result - - def tensorflow_layer_grad(op, grad): - """Thin wrapper for the gradient.""" - x = op.inputs[0] - return tensorflow_layer_grad_impl(x, grad) - - def tensorflow_layer(x): - """Implementation of the tensorflow call. - - Gradient in tensorflow is equivalent to the adjoint of the derivative - in ODL. - """ - x_shape = x.get_shape() - try: - n_x = int(x_shape[0]) - fixed_size = True - except TypeError: - n_x = x_shape[0] - fixed_size = False - - in_shape = (n_x,) + odl_op.domain.shape + (1,) - out_shape = (n_x,) + odl_op.range.shape + (1,) - - # Validate input shape - assert x_shape[1:] == odl_op.domain.shape + (1,) - - def _impl(x): - if fixed_size: - x_out_shape = out_shape - assert x.shape == in_shape - else: - x_out_shape = (x.shape[0],) + out_shape[1:] - assert x.shape[1:] == in_shape[1:] - - out = np.empty(x_out_shape, odl_op.range.dtype) - for i in range(x_out_shape[0]): - out[i] = np.asarray(odl_op(x[i, ..., 0]))[..., None] - return out - - if differentiable: - grad = tensorflow_layer_grad - else: - grad = None - - with ops.name_scope(name, "ODLTensorflowLayer", [x]) as name_call: - result = py_func(_impl, - [x], - [tf.float32], - name=name_call, - stateful=False, - grad=grad) - - # We must manually set the output shape since tensorflow cannot - # figure it out - result = result[0] - result.set_shape(out_shape) - return result - - return tensorflow_layer - - if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From b5634105054055208f9a9d2d1424bde9637bdbec Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 1 Aug 2017 14:07:46 +0200 Subject: [PATCH 57/70] MAINT: Bug fixes for product space --- odl/test/operator/tensor_ops_test.py | 20 ++++++++------------ odl/test/space/pspace_test.py | 2 +- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/odl/test/operator/tensor_ops_test.py b/odl/test/operator/tensor_ops_test.py index 478cacf12b1..15fd250eb7b 100644 --- a/odl/test/operator/tensor_ops_test.py +++ b/odl/test/operator/tensor_ops_test.py @@ -355,12 +355,11 @@ def test_pointwise_inner_adjoint(): testfunc = fspace.element(testarr) testfunc_pwinner_adj = pwinner.adjoint(testfunc) - assert all_almost_equal(testfunc_pwinner_adj, - true_inner_adj.reshape([1, -1])) + assert all_almost_equal(testfunc_pwinner_adj, true_inner_adj) out = vfspace.element() pwinner.adjoint(testfunc, out=out) - assert all_almost_equal(out, true_inner_adj.reshape([1, -1])) + assert all_almost_equal(out, true_inner_adj) # 3d fspace = odl.uniform_discr([0, 0], [1, 1], (2, 2), dtype=complex) @@ -380,12 +379,11 @@ def test_pointwise_inner_adjoint(): testfunc = fspace.element(testarr) testfunc_pwinner_adj = pwinner.adjoint(testfunc) - assert all_almost_equal(testfunc_pwinner_adj, - true_inner_adj.reshape([3, -1])) + assert all_almost_equal(testfunc_pwinner_adj, true_inner_adj) out = vfspace.element() pwinner.adjoint(testfunc, out=out) - assert all_almost_equal(out, true_inner_adj.reshape([3, -1])) + assert all_almost_equal(out, true_inner_adj) def test_pointwise_inner_adjoint_weighted(): @@ -407,12 +405,11 @@ def test_pointwise_inner_adjoint_weighted(): testfunc = fspace.element(testarr) testfunc_pwinner_adj = pwinner.adjoint(testfunc) - assert all_almost_equal(testfunc_pwinner_adj, - true_inner_adj.reshape([3, -1])) + assert all_almost_equal(testfunc_pwinner_adj, true_inner_adj) out = vfspace.element() pwinner.adjoint(testfunc, out=out) - assert all_almost_equal(out, true_inner_adj.reshape([3, -1])) + assert all_almost_equal(out, true_inner_adj) # Using different weighting in the inner product pwinner = PointwiseInner(vfspace, vecfield=array, weighting=[4, 8, 12]) @@ -424,12 +421,11 @@ def test_pointwise_inner_adjoint_weighted(): testfunc = fspace.element(testarr) testfunc_pwinner_adj = pwinner.adjoint(testfunc) - assert all_almost_equal(testfunc_pwinner_adj, - true_inner_adj.reshape([3, -1])) + assert all_almost_equal(testfunc_pwinner_adj, true_inner_adj) out = vfspace.element() pwinner.adjoint(testfunc, out=out) - assert all_almost_equal(out, true_inner_adj.reshape([3, -1])) + assert all_almost_equal(out, true_inner_adj) # ---- PointwiseSum ---- diff --git a/odl/test/space/pspace_test.py b/odl/test/space/pspace_test.py index c4b4e85cb41..bda68e6665c 100644 --- a/odl/test/space/pspace_test.py +++ b/odl/test/space/pspace_test.py @@ -38,7 +38,7 @@ def test_RxR(): # Check the basic properties assert len(HxH) == 2 - assert HxH.shape == (2,) + assert HxH.shape == (2, 2) assert HxH.size == 2 assert HxH.dtype == H.dtype assert HxH.spaces[0] is H From a25a575e4a89ec4a5e5bc5f19cb357993adcca4f Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 1 Aug 2017 15:04:33 +0200 Subject: [PATCH 58/70] ENH: Add TensorflowOperator --- odl/contrib/tensorflow/__init__.py | 3 + .../examples/tensorflow_operator.py | 142 ------------------ .../examples/tensorflow_operator_matrix.py | 53 +++++++ odl/contrib/tensorflow/operator.py | 109 ++++++++++++++ 4 files changed, 165 insertions(+), 142 deletions(-) delete mode 100644 odl/contrib/tensorflow/examples/tensorflow_operator.py create mode 100644 odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py create mode 100644 odl/contrib/tensorflow/operator.py diff --git a/odl/contrib/tensorflow/__init__.py b/odl/contrib/tensorflow/__init__.py index 5ae66c282a0..e81b62ef645 100644 --- a/odl/contrib/tensorflow/__init__.py +++ b/odl/contrib/tensorflow/__init__.py @@ -18,3 +18,6 @@ from .space import * __all__ += space.__all__ + +from .operator import * +__all__ += operator.__all__ diff --git a/odl/contrib/tensorflow/examples/tensorflow_operator.py b/odl/contrib/tensorflow/examples/tensorflow_operator.py deleted file mode 100644 index 5ca9e8aa855..00000000000 --- a/odl/contrib/tensorflow/examples/tensorflow_operator.py +++ /dev/null @@ -1,142 +0,0 @@ -"""Example of wrapping tensorflow to create an ODL operator.""" -import numpy as np -import matplotlib.pyplot as plt -from matplotlib.legend_handler import HandlerLine2D -import tensorflow as tf -import odl - - -sess = tf.InteractiveSession() - -def G(x): - if(np.array(x.shape).size==1): - return 1./(np.pi**.5)*np.exp(-x**2/2.) - elif(np.array(x.shape).size==2): - return 1./(np.pi**.5)*np.exp(-np.sum(x**2,axis=1)/2.) - else: - print('Error: x should be an MxN array') - return - -def G_sigma_m_x_m(x,sigma_m,x_m): - return G((x-x_m)/sigma_m) - -def fp(Arg,x): - # Arg = 4xM array - # Arg = [[alpha_j],[sigma_j],[x1_j],[x2_j]] - fp_temp = np.zeros(x.shape[0]) - M = Arg.shape[1] - alpha_j = Arg[0] - sigma_j = Arg[1] - x_j = Arg[2:].T - - for m in np.arange(M): - fp_temp = fp_temp + alpha_j[m]*G((x-x_j[m])/sigma_j[m]) - - return fp_temp - - -def R_G(p): - return tf.exp(-p ** 2. / 2.) - - -def R_G_sigma_m_x_m(theta, r, sigma_m, x_m, y_m): - return R_G((r - (tf.cos(theta) * x_m + tf.sin(theta) * y_m)) / sigma_m) - - -def R_fp(alpha, sigma, x, y, theta, r): - R_fp_temp = R_G_sigma_m_x_m(theta[:, None], r[:, None], - sigma[None, :], x[None, :], y[None, :]) - - R_fp_temp = alpha[None, :] * R_fp_temp - - return tf.reduce_sum(R_fp_temp, axis=-1) - - -class GaussianRayTransform(odl.Operator): - def __init__(self, domain, range, theta, p): - size = domain[0].size - self.alpha_ph = tf.placeholder(tf.float32, shape=[size]) - self.sigma_ph = tf.placeholder(tf.float32, shape=[size]) - self.x_ph = tf.placeholder(tf.float32, shape=[size]) - self.y_ph = tf.placeholder(tf.float32, shape=[size]) - theta = tf.constant(theta, dtype=tf.float32) - p = tf.constant(p, dtype=tf.float32) - - self.y = tf.placeholder(tf.float32, shape=range.shape) - - self.result = R_fp(self.alpha_ph, self.sigma_ph, self.x_ph, self.y_ph, theta, p) - self.gradient = tf.gradients(self.result, - [self.alpha_ph, self.sigma_ph, self.x_ph, self.y_ph], - [self.y]) - self.gradient = [gi * range.cell_volume for gi in self.gradient] - - odl.Operator.__init__(self, domain, range) - - def _call(self, x): - result = sess.run(self.result, - feed_dict={self.alpha_ph: np.asarray(x[0]), - self.sigma_ph: np.asarray(x[1]), - self.x_ph: np.asarray(x[2]), - self.y_ph: np.asarray(x[3])}) - - return result.reshape(self.range.shape).T - - def derivative(self, x): - op = self - - class GaussianRayTransformDerivative(odl.Operator): - @property - def adjoint(self): - class GaussianRayTransformDerivativeAdjoint(odl.Operator): - def _call(self, y): - result = sess.run(op.gradient, - feed_dict={op.alpha_ph: np.asarray(x[0]), - op.sigma_ph: np.asarray(x[1]), - op.x_ph: np.asarray(x[2]), - op.y_ph: np.asarray(x[3]), - op.y: np.asarray(y).T}) - - return result - - return GaussianRayTransformDerivativeAdjoint(self.range, - self.domain, - linear=True) - - return GaussianRayTransformDerivative(self.domain, self.range, - linear=True) - - -if __name__ == '__main__': - n = 201 - x,y = np.meshgrid(np.linspace(-2,2,n)*5,np.linspace(-2,2,n)*5) - X = np.row_stack((x.flatten(),y.flatten())).T - - M = 13 - r = np.linspace(0,2,M)*4 - x_theta = np.linspace(0,2.*np.pi,M) - x_ini = r*np.array([np.cos(x_theta),np.sin(x_theta)]) - #plt.plot(x_ini[0],x_ini[1],'.') - #plt.show() - alpha_ini = np.ones(M) - sigma_ini = np.ones(M) - Arg_ini = np.row_stack((alpha_ini,sigma_ini,x_ini)) - - theta,p = np.meshgrid(np.linspace(0,2.*np.pi,n),np.linspace(-20,20,n)) - Theta = theta.flatten() - P = p.flatten() - - import odl - rn = odl.rn(M) - domain = rn ** 4 - - ran = odl.uniform_discr([0, -20], [2 * np.pi, 20], [n, n]) - operator = GaussianRayTransform(domain, ran, Theta, P) - - x = domain.element([alpha_ini, sigma_ini, x_ini[0], x_ini[1]]) - - y = operator(x) - y.show() - - deriv = operator.derivative(x) - - adjy = deriv.adjoint(y) diff --git a/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py b/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py new file mode 100644 index 00000000000..3a1681fee6a --- /dev/null +++ b/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py @@ -0,0 +1,53 @@ +"""Example of how to convert a tensorflow layer to an ODL operator. + +In this example we take a tensorflow layer (given by matrix multiplication) +and convert it into an ODL operator. + +We also demonstrate that both the +""" + +import tensorflow as tf +import numpy as np +import odl +import odl.contrib.tensorflow + +sess = tf.InteractiveSession() + +matrix = np.array([[1, 2], + [0, 0], + [0, 1]], dtype='float32') + +# Define ODL operator +odl_op_pure = odl.MatrixOperator(matrix) + +# Define ODL operator using tensorflow +input_tensor = tf.placeholder(tf.float32, shape=[2]) +output_tensor = tf.reshape(tf.matmul(matrix, + tf.expand_dims(input_tensor, 1)), + [3]) + +odl_op_tensorflow = odl.contrib.tensorflow.TensorflowOperator(input_tensor, + output_tensor) + +# Define evaluation points +x = [1., 2.] +dx = [3., 4.] +dy = [1., 2., 3.] + +# Evaluate +print('Tensorflow eval: ', + odl_op_tensorflow(x)) +print('ODL eval: ', + odl_op_pure(x)) + +# Evaluate the derivative +print('Tensorflow derivative: ', + odl_op_tensorflow.derivative(x)(dx)) +print('ODL derivative: ', + odl_op_pure.derivative(x)(dx)) + +# Evaluate the adjoint of the derivative +print('Tensorflow adjoint of derivative: ', + odl_op_tensorflow.derivative(x).adjoint(dy)) +print('ODL adjoint of derivative: ', + odl_op_pure.derivative(x).adjoint(dy)) diff --git a/odl/contrib/tensorflow/operator.py b/odl/contrib/tensorflow/operator.py new file mode 100644 index 00000000000..ffb627c7c2a --- /dev/null +++ b/odl/contrib/tensorflow/operator.py @@ -0,0 +1,109 @@ +# Copyright 2014-2017 The ODL development group +# +# This file is part of ODL. +# +# ODL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ODL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with ODL. If not, see . + +"""Utilities for converting ODL spaces to tensorflow layers.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import tensorflow as tf +import numpy as np +import odl + + +__all__ = ('TensorflowOperator',) + + +class TensorflowOperator(odl.Operator): + def __init__(self, input_tensor, output_tensor, linear=False, sess=None): + self.input_tensor = input_tensor + self.output_tensor = output_tensor + + # TODO: Fix with tensors + domain = odl.fn(np.prod(input_tensor.shape.as_list()), + dtype=input_tensor.dtype.as_numpy_dtype) + range = odl.fn(np.prod(output_tensor.shape.as_list()), + dtype=output_tensor.dtype.as_numpy_dtype) + + self.dx = tf.placeholder(input_tensor.dtype, + shape=input_tensor.shape) + self.dy = tf.placeholder(output_tensor.dtype, + shape=output_tensor.shape) + + adjoint_of_derivative_tensor = tf.gradients( + self.output_tensor, [self.input_tensor], [self.dy])[0] + self.adjoint_of_derivative_tensor = (range.weighting.const * + adjoint_of_derivative_tensor) + + # Since tensorflow does not support forward differentiation, use trick + # that adjoint of the derivative of adjoint of the derivative is simply + # the derivative. + derivative_tensor = tf.gradients( + adjoint_of_derivative_tensor, [self.dy], [self.dx])[0] + self.derivative_tensor = (range.weighting.const * + derivative_tensor) + + if sess is None: + self.sess = tf.get_default_session() + else: + self.sess = sess + + odl.Operator.__init__(self, domain, range, linear=linear) + + def _call(self, x): + result = self.sess.run(self.output_tensor, + feed_dict={self.input_tensor: np.asarray(x)}) + + return result + + def derivative(self, x): + op = self + + class TensorflowOperatorDerivative(odl.Operator): + def _call(self, dx): + result = op.sess.run(op.derivative_tensor, + feed_dict={op.input_tensor: np.asarray(x), + op.dx: np.asarray(dx)}) + + return result + + @property + def adjoint(self): + class TensorflowOperatorDerivativeAdjoint(odl.Operator): + def _call(self, y): + result = op.sess.run( + op.adjoint_of_derivative_tensor, + feed_dict={op.input_tensor: np.asarray(x), + op.dy: np.asarray(y)}) + + return result + + return TensorflowOperatorDerivativeAdjoint(self.range, + self.domain, + linear=True) + + return TensorflowOperatorDerivative(self.domain, + self.range, + linear=True) + + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() From 880922524332e8cab8ff22aef7ef5c0cc2e5ade2 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Thu, 3 Aug 2017 15:46:07 +0200 Subject: [PATCH 59/70] MAINT: Improve tensorflow examples --- .../examples/tensorflow_layer_matrix.py | 25 +++++++++---- .../examples/tensorflow_layer_productspace.py | 12 +++++-- .../tensorflow_layer_ray_transform.py | 6 +++- .../examples/tensorflow_operator_matrix.py | 2 -- .../examples/tensorflow_tomography.py | 36 +++++++++++++------ odl/contrib/tensorflow/layer.py | 6 ++-- 6 files changed, 63 insertions(+), 24 deletions(-) diff --git a/odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py b/odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py index 41f32998f13..e40b30bcd86 100644 --- a/odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py +++ b/odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py @@ -1,8 +1,17 @@ -"""Example of how to convert an ODL operator to a tensorflow layer.""" +"""Example of how to convert an ODL operator to a tensorflow layer. + +In this example we take an ODL operator given by a `MatrixOperator` and +convert it into a tensorflow layer that can be used inside any tensorflow +computational graph. + +We also demonstrate that we can compute the "gradients" properly using the +adjoint of the derivative. +""" import tensorflow as tf import numpy as np import odl +import odl.contrib.tensorflow sess = tf.InteractiveSession() @@ -22,17 +31,21 @@ # Create tensorflow layer from odl operator odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( - odl_op, 'MatrixOperator') + odl_op, 'MatrixOperator') y_tf = odl_op_layer(x_tf) # Evaluate using tensorflow -print(y_tf.eval().ravel()) +print('Tensorflow eval: ', + y_tf.eval().ravel()) # Compare result with pure ODL -print(odl_op(x)) +print('ODL eval: ', + odl_op(x)) # Evaluate the adjoint of the derivative, called gradient in tensorflow -print(tf.gradients(y_tf, [x_tf], z_tf)[0].eval().ravel()) +print('Tensorflow adjoint of derivative: ', + tf.gradients(y_tf, [x_tf], z_tf)[0].eval().ravel()) # Compare result with pure ODL -print(odl_op.derivative(x).adjoint(z)) +print('ODL adjoint of derivative: ', + odl_op.derivative(x).adjoint(z)) diff --git a/odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py b/odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py index 7664ed619bd..00f4ee1e51b 100644 --- a/odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py +++ b/odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py @@ -1,8 +1,12 @@ -"""Example of how to convert an ODL operator to a tensorflow layer.""" +"""Example of how to convert an ODL operator to a tensorflow layer. + +This example is similar to ``tensorflow_layer_matrix``, but demonstrates how +to handle product-spaces. +""" import tensorflow as tf -import numpy as np import odl +import odl.contrib.tensorflow sess = tf.InteractiveSession() @@ -25,14 +29,18 @@ y_tf = odl_op_layer(x_tf) # Evaluate using tensorflow +print('Tensorflow eval:') print(y_tf.eval().squeeze()) # Compare result with pure ODL +print('ODL eval') print(odl_op(x)) # Evaluate the adjoint of the derivative, called gradient in tensorflow scale = 1 / space.cell_volume +print('Tensorflow adjoint of derivative (gradients):') print(tf.gradients(y_tf, [x_tf], z_tf)[0].eval().squeeze() * scale) # Compare result with pure ODL +print('ODL adjoint of derivative:') print(odl_op.derivative(x).adjoint(z)) diff --git a/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py b/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py index 3197e214200..0ccc3ef4bd2 100644 --- a/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py +++ b/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py @@ -1,4 +1,8 @@ -"""Example of how to convert an ODL ray transform to a tensorflow layer.""" +"""Example of how to convert a RayTransform operator to a tensorflow layer. + +This example is similar to ``tensorflow_layer_matrix``, but demonstrates how +more advanced operators, such as a ray transform, can be handled. +""" import tensorflow as tf import numpy as np diff --git a/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py b/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py index 3a1681fee6a..e90c81c738a 100644 --- a/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py +++ b/odl/contrib/tensorflow/examples/tensorflow_operator_matrix.py @@ -2,8 +2,6 @@ In this example we take a tensorflow layer (given by matrix multiplication) and convert it into an ODL operator. - -We also demonstrate that both the """ import tensorflow as tf diff --git a/odl/contrib/tensorflow/examples/tensorflow_tomography.py b/odl/contrib/tensorflow/examples/tensorflow_tomography.py index 23a377c2d14..242e008e589 100644 --- a/odl/contrib/tensorflow/examples/tensorflow_tomography.py +++ b/odl/contrib/tensorflow/examples/tensorflow_tomography.py @@ -1,9 +1,16 @@ -"""Example of how to solve a simple tomography problem using tensorflow.""" +"""Example of how to solve a simple tomography problem using tensorflow. + +In this example, we solve the TV regularized tomography problem:: + + min_x ||A(x) - b||_2^2 + 50 * ||grad(x)||_1 + +using a gradient descent method (ADAM). +""" import tensorflow as tf import numpy as np import odl - +import odl.contrib.tensorflow sess = tf.InteractiveSession() @@ -12,25 +19,32 @@ dtype='float32') geometry = odl.tomo.parallel_beam_geometry(space) ray_transform = odl.tomo.RayTransform(space, geometry) +grad = odl.Gradient(space) + +# Create data phantom = odl.phantom.shepp_logan(space, True) data = ray_transform(phantom) +noisy_data = data + odl.phantom.white_noise(data.space) -# Create tensorflow layer from odl operator -odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( - ray_transform, 'RayTransform') -x = tf.Variable(tf.constant(0.0, shape=space.shape), name="x") +# Create tensorflow layers from odl operators +ray_transform_layer = odl.contrib.tensorflow.as_tensorflow_layer( + ray_transform, name='RayTransform') +grad_layer = odl.contrib.tensorflow.as_tensorflow_layer( + grad, name='Gradient') +x = tf.Variable(tf.zeros(shape=space.shape), name="x") # Create constant right hand side -y = tf.constant(np.asarray(data)) +y = tf.constant(np.asarray(noisy_data)) # Add empty axes for batch and channel x = x[None, ..., None] y = y[None, ..., None] -# Define l2 loss function -loss = tf.nn.l2_loss(odl_op_layer(x) - y) +# Define loss function +loss = (tf.reduce_sum((ray_transform_layer(x) - y) ** 2) + + 50 * tf.reduce_sum(tf.abs(grad_layer(x)))) -# Train using the adam optimizer +# Train using the ADAM optimizer optimizer = tf.train.AdamOptimizer(1e-1).minimize(loss) # Initialize all TF variables @@ -39,7 +53,7 @@ # Solve with an ODL callback to see what happens callback = odl.solvers.CallbackShow() -for i in range(100): +for i in range(200): sess.run(optimizer) callback(space.element(x.eval())) diff --git a/odl/contrib/tensorflow/layer.py b/odl/contrib/tensorflow/layer.py index 4b49a576071..822d03a533e 100644 --- a/odl/contrib/tensorflow/layer.py +++ b/odl/contrib/tensorflow/layer.py @@ -32,7 +32,7 @@ __all__ = ('as_tensorflow_layer',) -def as_tensorflow_layer(odl_op, default_name='ODLOperator', +def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): """Convert `Operator` or `Functional` into tensorflow layer. @@ -40,7 +40,7 @@ def as_tensorflow_layer(odl_op, default_name='ODLOperator', ---------- odl_op : `Operator` or `Functional` The operator that should be wrapped to a tensorflow layer. - default_name : str + name : str Default name for tensorflow layers created. differentiable : boolean True if the layer should be differentiable, in which case ``odl_op`` @@ -64,6 +64,8 @@ def as_tensorflow_layer(odl_op, default_name='ODLOperator', The `dtype` of the tensor is the same as the respective ODL spaces. """ + default_name = name + def py_func(func, inp, Tout, stateful=True, name=None, grad=None): """Define custom py_func which takes also a grad op as argument. From e6db7e041bf8ef4bba75cd3acc9937f22fa517d8 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 18 Aug 2017 14:17:28 +0200 Subject: [PATCH 60/70] DOC: Remove overview_functionality stub --- doc/source/guide/overview_functionality.rst | 103 -------------------- 1 file changed, 103 deletions(-) delete mode 100644 doc/source/guide/overview_functionality.rst diff --git a/doc/source/guide/overview_functionality.rst b/doc/source/guide/overview_functionality.rst deleted file mode 100644 index 25071a0aa0f..00000000000 --- a/doc/source/guide/overview_functionality.rst +++ /dev/null @@ -1,103 +0,0 @@ -.. _overview_functionality: - -######################### -Overview of functionality -######################### - - -This section gives a short overview of the functionality available in ODL. - - -What is vectorization? -====================== - -In general, :term:`vectorization` means that a function can be evaluated on a whole array of values -at once instead of looping over individual entries. This is very important for performance in an -interpreted language like python, since loops are usually very slow compared to compiled languages. - -Technically, vectorization in Numpy works through the `Universal functions (ufunc)`_ interface. It -is fast because all loops over data are implemented in C, and the resulting implementations are -exposed to Python for each function individually. - - -How to use Numpy's ufuncs? -========================== - -The easiest way to write fast custom mathematical functions in Python is to use the -`available ufuncs`_ and compose them to a new function:: - - def gaussian(x): - # Negation, powers and scaling are vectorized, of course. - return np.exp(-x ** 2 / 2) - - def step(x): - # np.where checks the condition in the first argument and - # returns the second for `True`, otherwise the third. The - # last two arguments can be arrays, too. - # Note that also the comparison operation is vectorized. - return np.where(x[0] <= 0, 0, 1) - -This should cover a very large range of useful functions already (basic arithmetic is vectorized, -too!). An even larger list of `special functions`_ are available in the Scipy package. - - -Usage in ODL -============ - -Python functions are in most cases used as input to a discretization process. For example, we may -want to discretize a two-dimensional Gaussian function: - - >>> def gaussian2(x): - ... return np.exp(-(x[0] ** 2 + x[1] ** 2) / 2) - -on the rectangle [-5, 5] x [-5, 5] with 100 pixels in each -dimension. The code for this is simply - - >>> # Note that the minimum and maxiumum coordinates are given as - >>> # vectors, not one interval at a time. - >>> discr = odl.uniform_discr([-5, -5], [5, 5], (100, 100)) - - >>> # This creates an element in the discretized space ``discr`` - >>> gaussian_discr = discr.element(gaussian2) - -What happens behind the scenes is that ``discr`` creates a :term:`discretization` object which -has a built-in method ``element`` to turn continuous functions into discrete arrays by evaluating -them at a set of grid points. In the example above, this grid is a uniform sampling of the rectangle -by 100 points per dimension. - -To make this process fast, ``element`` assumes that the function is written in a way that not only -supports vectorization, but also guarantees that the output has the correct shape. The function -receives a :term:`meshgrid` tuple as input, in the above case consisting of two vectors:: - - >>> mesh = discr.meshgrid - >>> mesh[0].shape - (100, 1) - >>> mesh[1].shape - (1, 100) - -When inserted into the function, the final shape of the output is determined by Numpy's -`broadcasting rules`_. For the Gaussian function, Numpy will conclude that the output shape must -be ``(100, 100)`` since the arrays in ``mesh`` are added after squaring. This size is the same -as expected by the discretization. - -If a function does not use all components of the input, ODL tries to broadcast the result to the shape of the discretized space:: - - >>> def gaussian_const_x0(x): - ... return np.exp(-x[1] ** 2 / 2) # no x[0] -> broadcasting - - >>> gaussian_const_x0(mesh).shape - (1, 100) - >>> discr.element(gaussian_const_x0).shape - (100, 100) - - -Further reading -=============== - -`Scipy Lecture notes on Numpy `_ - - -.. _Universal functions (ufunc): http://docs.scipy.org/doc/numpy/reference/ufuncs.html -.. _available ufuncs: http://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs -.. _special functions: http://docs.scipy.org/doc/scipy/reference/special.html -.. _broadcasting rules: http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html From c397833bac68a2fda152b0f9d4d82928c72f9a7d Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 18 Aug 2017 14:55:48 +0200 Subject: [PATCH 61/70] MAINT: Doc improvements to adam solver --- odl/solvers/smooth/gradient.py | 36 +++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/odl/solvers/smooth/gradient.py b/odl/solvers/smooth/gradient.py index aebbf065f67..e6203650a76 100644 --- a/odl/solvers/smooth/gradient.py +++ b/odl/solvers/smooth/gradient.py @@ -118,41 +118,45 @@ def adam(f, x, learning_rate=1e-3, beta1=0.9, beta2=0.999, eps=1e-8, .. math:: \min f(x) - The algorithm is intended for unconstrained problems. + where :math:`f` is a differentiable functional. - The algorithm is described in - `Adam: A Method for Stochastic Optimization - `_. All parameter names are taken from - that article. + The algorithm is described in [KB2015] (`arxiv + `_). All parameter names and default + valuesare taken from the article. Parameters ---------- f : `Functional` Goal functional. Needs to have ``f.gradient``. x : ``f.domain`` element - Starting point of the iteration - learning_rate : float, optional + Starting point of the iteration, updated in place. + learning_rate : positive float, optional Step length of the method. - beta1 : float, optional + beta1 : float in [0, 1), optional Update rate for first order moment estimate. - beta2 : float, optional + beta2 : float in [0, 1), optional Update rate for second order moment estimate. - eps : float, optional + eps : positive float, optional A small constant for numerical stability. maxiter : int, optional Maximum number of iterations. - tol : float, optional + tol : positive float, optional Tolerance that should be used for terminating the iteration. callback : callable, optional - Object executing code per iteration, e.g. plotting each iterate + Object executing code per iteration, e.g. plotting each iterate. See Also -------- - odl.solvers.smooth.gradient.steepest_descent : simple steepest descent + odl.solvers.smooth.gradient.steepest_descent : Simple gradient descent. odl.solvers.iterative.iterative.landweber : - Optimized solver for the case ``f(x) = ||Ax - b||_2^2`` + Optimized solver for the case ``f(x) = ||Ax - b||_2^2``. odl.solvers.iterative.iterative.conjugate_gradient : - Optimized solver for the case ``f(x) = x^T Ax - 2 x^T b`` + Optimized solver for the case ``f(x) = x^T Ax - 2 x^T b``. + + References + ---------- + [KB2015] Kingma, D P and Ba, J. + *Adam: A Method for Stochastic Optimization*, ICLR 2015. """ grad = f.gradient if x not in grad.domain: @@ -174,7 +178,7 @@ def adam(f, x, learning_rate=1e-3, beta1=0.9, beta2=0.999, eps=1e-8, step = learning_rate * np.sqrt(1 - beta2) / (1 - beta1) - x.lincomb(1, x, -step, m / np.sqrt(v + eps)) + x.lincomb(1, x, -step, m / (np.sqrt(v) + eps)) if callback is not None: callback(x) From 6d064884cd7556ad314b000cc752a5a84b0a9d76 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Fri, 18 Aug 2017 15:06:22 +0200 Subject: [PATCH 62/70] MAINT: Improve doc of ProductSpace.element.__array__ --- odl/space/pspace.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/odl/space/pspace.py b/odl/space/pspace.py index 4dc5d7a7aed..86ab7b3b3ed 100644 --- a/odl/space/pspace.py +++ b/odl/space/pspace.py @@ -633,7 +633,7 @@ def parts(self): def shape(self): """Number of spaces per axis.""" return self.space.shape - + @property def size(self): """Number of factors of this element's space.""" @@ -715,15 +715,20 @@ def __setitem__(self, indices, values): def __array__(self): """An array representation of ``self``. - + Only available if `is_power_space` is True. - + + The ordering is such that it commutes with indexing:: + + np.asarray(self[ind]) == np.asarray(self)[ind] + Examples -------- - >>> spc = odl.ProductSpace(odl.rn(2), 2) + >>> spc = odl.ProductSpace(odl.rn(2), 3) >>> x = spc.one() >>> np.asarray(x) array([[ 1., 1.], + [ 1., 1.], [ 1., 1.]]) """ if not self.space.is_power_space: From c7e0ca5fbf6c4956e8b6cb30c4d435a2be716740 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:11:45 +0200 Subject: [PATCH 63/70] ENH: Add zscore --- odl/util/numerics.py | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/odl/util/numerics.py b/odl/util/numerics.py index 3f5f9a33794..052f63f55fd 100644 --- a/odl/util/numerics.py +++ b/odl/util/numerics.py @@ -18,7 +18,8 @@ from odl.util.normalize import normalized_scalar_param_list, safe_int_conv -__all__ = ('apply_on_boundary', 'fast_1d_tensor_mult', 'resize_array') +__all__ = ('apply_on_boundary', 'fast_1d_tensor_mult', 'resize_array', + 'zscore') _SUPPORTED_RESIZE_PAD_MODES = ('constant', 'symmetric', 'periodic', @@ -802,6 +803,42 @@ def _apply_padding(lhs_arr, rhs_arr, offset, pad_mode, direction): working_slc[axis] = intersec_slc[axis] +def zscore(arr): + """Return arr normalized with mean 0 and unit variance. + + If the input has 0 variance, the result will also have 0 variance. + + Parameters + ---------- + arr : array-like + + Returns + ------- + zscore : array-like + + Examples + -------- + Compute the z score for a small array: + + >>> result = zscore([1, 0]) + >>> result + array([ 1., -1.]) + >>> np.mean(result) + 0.0 + >>> np.std(result) + 1.0 + + Does not re-scale in case the input is constant (has 0 variance): + + >>> zscore([1, 1]) + array([ 0., 0.]) + """ + arr = arr - np.mean(arr) + std = np.std(arr) + if std != 0: + arr /= std + return arr + if __name__ == '__main__': # pylint: disable=wrong-import-position from odl.util.testutils import run_doctests From 38e3ab6b04af4a1c3dd78c7a398b672d62afacac Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:12:43 +0200 Subject: [PATCH 64/70] ENH: Add util/statistics to contrib/fom --- odl/contrib/fom/__init__.py | 3 + odl/contrib/fom/supervised.py | 61 ++++++++- odl/contrib/fom/test/test_unsupervised.py | 57 +++++++++ odl/contrib/fom/unsupervised.py | 65 ++++++++++ odl/util/__init__.py | 3 - odl/util/statistics.py | 147 ---------------------- 6 files changed, 185 insertions(+), 151 deletions(-) create mode 100644 odl/contrib/fom/test/test_unsupervised.py create mode 100644 odl/contrib/fom/unsupervised.py delete mode 100644 odl/util/statistics.py diff --git a/odl/contrib/fom/__init__.py b/odl/contrib/fom/__init__.py index 9dffd8bfba9..f350af714a6 100644 --- a/odl/contrib/fom/__init__.py +++ b/odl/contrib/fom/__init__.py @@ -12,3 +12,6 @@ from .supervised import * __all__ += supervised.__all__ + +from .unsupervised import * +__all__ += unsupervised.__all__ diff --git a/odl/contrib/fom/supervised.py b/odl/contrib/fom/supervised.py index 2710c66d042..7b246dbeb3a 100644 --- a/odl/contrib/fom/supervised.py +++ b/odl/contrib/fom/supervised.py @@ -16,7 +16,7 @@ __all__ = ('mean_squared_error', 'mean_absolute_error', 'mean_value_difference', 'standard_deviation_difference', - 'range_difference', 'blurring', 'false_structures', 'ssim') + 'range_difference', 'blurring', 'false_structures', 'ssim', 'psnr') def mean_squared_error(data, ground_truth, mask=None, normalized=False): @@ -573,3 +573,62 @@ def smoothen(img): return 0.5 - ssim / 2 else: return ssim + + +def psnr(data, ground_truth, normalize=False): + """Return the Peak Signal-to-Noise Ratio. + + Parameters + ---------- + data : `FnBaseVector` + Input data or reconstruction. + ground_truth : `FnBaseVector` + Reference to compare ``data`` to. + normalize : bool + If true, normalizes ``data`` and ``ground_truth`` to have the same mean + and variance before comparison. + + Returns + ------- + psnr : float + + Examples + -------- + Compute the PSNR for two vectors: + + >>> spc = odl.rn(5) + >>> data = spc.element([1, 1, 1, 1, 1]) + >>> ground_truth = spc.element([1, 1, 1, 1, 2]) + >>> result = psnr(data, ground_truth) + >>> print('{:.3f}'.format(result)) + 6.021 + + If data == ground_truth, the result is positive infinity: + + >>> psnr(ground_truth, ground_truth) + inf + + With ``normalize=True``, internal scaling and constant offsets are ignored: + + >>> (psnr(data, ground_truth, normalize=True) == + ... psnr(data, 3 + 4 * ground_truth, normalize=True)) + True + """ + if normalize: + data = odl.util.zscore(data) + ground_truth = odl.util.zscore(ground_truth) + + mse_result = mean_squared_error(data, ground_truth) + max_true = np.max(np.abs(ground_truth)) + + if mse_result == 0: + return np.inf + elif max_true == 0: + return -np.inf + else: + return 20 * np.log10(max_true) - 10 * np.log10(mse_result) + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/contrib/fom/test/test_unsupervised.py b/odl/contrib/fom/test/test_unsupervised.py new file mode 100644 index 00000000000..c713430cff8 --- /dev/null +++ b/odl/contrib/fom/test/test_unsupervised.py @@ -0,0 +1,57 @@ +# Copyright 2014-2017 The ODL contributors +# +# This file is part of ODL. +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. + +"""Tests for unsupervised FoMs.""" + +import numpy as np +import pytest +import odl +import odl.contrib.fom + + +def test_estimate_noise_std_constant_1d(): + """Verify ``estimate_noise_std(0) == 0``""" + img = np.zeros(10) + result = odl.contrib.fom.estimate_noise_std(img) + assert pytest.approx(result) == 0.0 + + +def test_estimate_noise_std_normal_1d(): + """Verify ``estimate_noise_std(N(0, 1)) == 1`` in 1d.""" + img = np.random.randn(10) + result = odl.contrib.fom.estimate_noise_std(img) + expected = np.std(img) + assert pytest.approx(result, abs=0.2) == expected + + +def test_estimate_noise_std_normal_2d(): + """Verify ``estimate_noise_std(N(0, 1)) == 1`` in 2d.""" + img = np.random.randn(10, 10) + result = odl.contrib.fom.estimate_noise_std(img) + expected = np.std(img) + assert pytest.approx(result, abs=0.2) == expected + + +def test_estimate_noise_std_normal_4d(): + """Verify ``estimate_noise_std(N(0, 1)) == 1`` in 4d.""" + img = np.random.randn(5, 5, 5, 5) + result = odl.contrib.fom.estimate_noise_std(img) + expected = np.std(img) + assert pytest.approx(result, abs=0.2) == expected + + +def test_estimate_noise_std_normal_large_1d(): + """Verify ``estimate_noise_std(N(0, 1)) == 1`` with low error.""" + img = np.random.randn(100000) + result = odl.contrib.fom.estimate_noise_std(img) + expected = np.std(img) + assert pytest.approx(result, abs=0.01) == expected + + +if __name__ == '__main__': + pytest.main([str(__file__.replace('\\', '/')), '-v']) diff --git a/odl/contrib/fom/unsupervised.py b/odl/contrib/fom/unsupervised.py new file mode 100644 index 00000000000..b9256d40627 --- /dev/null +++ b/odl/contrib/fom/unsupervised.py @@ -0,0 +1,65 @@ +# Copyright 2014-2017 The ODL contributors +# +# This file is part of ODL. +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. + +"""Figures of Merit (FOMs) for measuring image quality without a reference.""" + +import numpy as np + +__all__ = ('estimate_noise_std',) + + +def estimate_noise_std(img): + """Estimate standard deviation of noise in ``img``. + + The algorithm, given in [Immerkaer1996], estimates the noise in an image + by + + Parameters + ---------- + img : array-like + + Returns + ------- + noise : float + + Examples + -------- + Create image with noise 1.0, verify result + + >>> img = np.random.randn(10, 10) + >>> result = estimate_noise_std(img) # should be about 1 + + Also works with higher dimensional arrays + + >>> img = np.random.randn(3, 3, 3) + >>> result = estimate_noise_std(img) + + References + ---------- + [Immerkaer1996] Immerkaer, J. *Fast Noise Variance Estimation*. + Computer Vision and Image Understanding, 1996. + """ + import scipy.signal + import functools + img = np.asarray(img, dtype='float') + + M = functools.reduce(np.add.outer, [[-1, 2, -1]] * img.ndim) + + convolved = scipy.signal.fftconvolve(img, M, mode='valid') + conv_var = np.sum(convolved ** 2) + + scale = np.sum(np.square(M)) * convolved.size + sigma = np.sqrt(conv_var / scale) + + return sigma + + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/util/__init__.py b/odl/util/__init__.py index 54fff755bfd..9f57e6dcd00 100644 --- a/odl/util/__init__.py +++ b/odl/util/__init__.py @@ -27,9 +27,6 @@ from .numerics import * __all__ += numerics.__all__ -from .statistics import * -__all__ += statistics.__all__ - from .vectorization import * __all__ += vectorization.__all__ diff --git a/odl/util/statistics.py b/odl/util/statistics.py deleted file mode 100644 index f9d15726021..00000000000 --- a/odl/util/statistics.py +++ /dev/null @@ -1,147 +0,0 @@ -# Copyright 2014-2016 The ODL development group -# -# This file is part of ODL. -# -# ODL is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ODL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODL. If not, see . - -"""Utilities for computing statistics on images.""" - -# Imports for common Python 2/3 codebase -from __future__ import print_function, division, absolute_import -from future import standard_library -standard_library.install_aliases() - -import numpy as np - - -__all__ = ('mse', 'psnr', 'estimate_noise') - - -def mse(true, noisy): - """Return the Mean Squared Errror. - - Parameters - ---------- - true : array-like - noisy : array-like - - Returns - ------- - mse : float - - Examples - -------- - >>> true = [1, 1, 1, 1, 1] - >>> noisy = [1, 1, 1, 1, 2] - >>> result = mse(true, noisy) - >>> print('{:.3f}'.format(result)) - 0.200 - """ - true = np.asarray(true) - noisy = np.asarray(noisy) - return np.mean((true - noisy) ** 2) - - -def psnr(true, noisy, normalize=False): - """Return the Peak signal-to-noise ratio. - - Parameters - ---------- - true : array-like - noisy : array-like - normalize : bool - If true, normalizes the true and noisy signal to have the same mean - and variance before comparison. - - Returns - ------- - psnr : float - - Examples - -------- - >>> true = [1, 1, 1, 1, 1] - >>> noisy = [1, 1, 1, 1, 2] - >>> result = psnr(true, noisy) - >>> print('{:.3f}'.format(result)) - 6.990 - - If true == noisy, the result is infinite - - >>> psnr([1, 1], [1, 1]) - inf - - If `true == 0` but `noisy != 0`, the result is negative infinity - - >>> psnr(0, 1) - -inf - """ - true = np.asarray(true) - noisy = np.asarray(noisy) - - if normalize: - noisy = noisy + np.mean(true - noisy) - noisy = noisy * np.std(true) / np.std(noisy) - - mse_result = mse(true, noisy) - max_true = np.max(np.abs(true)) - - if mse_result == 0: - return np.inf - elif max_true == 0: - return -np.inf - else: - return 20 * np.log10(max_true) - 10 * np.log10(mse_result) - - -def estimate_noise(img): - """Estimate noise level in image. - - Parameters - ---------- - img : array-lik - - Returns - ------- - noise : float - estimate of standard deviation in image - - Examples - -------- - Create image with noise 1.0, verify result - - >>> img = np.random.randn(10, 10) - >>> noise = estimate_noise(img) # should be about 1 - - See Also - -------- - Reference: J. Immerkaer, “Fast Noise Variance Estimation” - """ - import scipy.signal - import functools - img = np.asarray(img, dtype='float') - - M = functools.reduce(np.add.outer, [[-1, 2, -1]] * img.ndim) - - convolved = scipy.signal.convolve(img, M, 'valid') - conv_var = np.sum(np.square(convolved)) - - sigma = np.sqrt(conv_var / (np.sum(M**2) * np.prod(convolved.shape))) - - return sigma - - -if __name__ == '__main__': - # pylint: disable=wrong-import-position - from odl.util.testutils import run_doctests - run_doctests() From 0acdc3cb2419055491dd6dd660e5634cec0b5b88 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:22:46 +0200 Subject: [PATCH 65/70] MAINT: Update licence of contrib/tensorflow --- odl/contrib/tensorflow/layer.py | 17 ++++------------- odl/contrib/tensorflow/operator.py | 17 ++++------------- odl/contrib/tensorflow/space.py | 17 ++++------------- 3 files changed, 12 insertions(+), 39 deletions(-) diff --git a/odl/contrib/tensorflow/layer.py b/odl/contrib/tensorflow/layer.py index 822d03a533e..4642689efb1 100644 --- a/odl/contrib/tensorflow/layer.py +++ b/odl/contrib/tensorflow/layer.py @@ -1,19 +1,10 @@ -# Copyright 2014-2017 The ODL development group +# Copyright 2014-2017 The ODL contributors # # This file is part of ODL. # -# ODL is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ODL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODL. If not, see . +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. """Utilities for converting ODL operators to tensorflow layers.""" diff --git a/odl/contrib/tensorflow/operator.py b/odl/contrib/tensorflow/operator.py index ffb627c7c2a..af247426c0c 100644 --- a/odl/contrib/tensorflow/operator.py +++ b/odl/contrib/tensorflow/operator.py @@ -1,19 +1,10 @@ -# Copyright 2014-2017 The ODL development group +# Copyright 2014-2017 The ODL contributors # # This file is part of ODL. # -# ODL is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ODL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODL. If not, see . +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. """Utilities for converting ODL spaces to tensorflow layers.""" diff --git a/odl/contrib/tensorflow/space.py b/odl/contrib/tensorflow/space.py index 25bb9daf217..e8811ef79ee 100644 --- a/odl/contrib/tensorflow/space.py +++ b/odl/contrib/tensorflow/space.py @@ -1,19 +1,10 @@ -# Copyright 2014-2017 The ODL development group +# Copyright 2014-2017 The ODL contributors # # This file is part of ODL. # -# ODL is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ODL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODL. If not, see . +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. """Utilities for converting ODL spaces to tensorflow layers.""" From 1737e0e0509f4f74c879907aa18e2b147f537a1f Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:22:56 +0200 Subject: [PATCH 66/70] DOC: Add readme to contrib/tensorflow --- odl/contrib/tensorflow/README.md | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 odl/contrib/tensorflow/README.md diff --git a/odl/contrib/tensorflow/README.md b/odl/contrib/tensorflow/README.md new file mode 100644 index 00000000000..f826a4cea7b --- /dev/null +++ b/odl/contrib/tensorflow/README.md @@ -0,0 +1,22 @@ +# Tensorflow + +This package contains ODL functionality related to [TensorFlow](https://www.tensorflow.org/). + +## Content + +* `TensorflowOperator` in [operator.py](operator.py) wraps a tensorflow expression into an ODL operator. + This allows using tensorflow neural networks as operators in ODL. +* `as_tensorflow_layer` in [layer.py](layer.py) wraps an ODL operator into a tensorflow layer. + This allows using arbitrary ODL operators inside tensorflow neural networks. +* `TensorflowSpace` in [space.py](space.py) is a `FnBase`-type space which uses tensorflow as a backend. + +## Example usage + +The [examples](examples) folder contains example on how to use the above functionality. +Specifically + +* [tensorflow_layer_matrix.py](examples/tensorflow_layer_matrix.py) shows how an ODL `MatrixOperator` can be converted to a tensorflow layer. +* [tensorflow_layer_productspace.py](examples/tensorflow_layer_productspace.py) shows how an ODL operator acting on `ProductSpace`s can be converted to a tensorflow layer. +* [tensorflow_layer_ray_transform.py](examples/tensorflow_layer_ray_transform.py) shows how a `RayTransform` can be converted to a tensorflow layer. +* [tensorflow_operator_matrix.py](examples/tensorflow_operator_matrix.py) shows how `tf.matmul` can be used as a ODL operator. +* [tensorflow_tomography.py](examples/tensorflow_tomography.py) shows how tensorflow optimizers can be used with ODL operators to solve inverse problems. From 127814efe81621b9a0cee5453f00b277c10069c3 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:41:52 +0200 Subject: [PATCH 67/70] MAINT: imporovements to contrib/tensorflow --- .../tensorflow_layer_ray_transform.py | 1 + odl/contrib/tensorflow/layer.py | 111 ++++++++++-------- 2 files changed, 61 insertions(+), 51 deletions(-) diff --git a/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py b/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py index 0ccc3ef4bd2..17a7d772cce 100644 --- a/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py +++ b/odl/contrib/tensorflow/examples/tensorflow_layer_ray_transform.py @@ -7,6 +7,7 @@ import tensorflow as tf import numpy as np import odl +import odl.contrib.tensorflow sess = tf.InteractiveSession() diff --git a/odl/contrib/tensorflow/layer.py b/odl/contrib/tensorflow/layer.py index 4642689efb1..361c2b30004 100644 --- a/odl/contrib/tensorflow/layer.py +++ b/odl/contrib/tensorflow/layer.py @@ -25,43 +25,48 @@ def as_tensorflow_layer(odl_op, name='ODLOperator', differentiable=True): - """Convert `Operator` or `Functional` into tensorflow layer. + """Convert `Operator` or `Functional` to a tensorflow layer. Parameters ---------- odl_op : `Operator` or `Functional` - The operator that should be wrapped to a tensorflow layer. + The operator that should be wrapped as a tensorflow layer. name : str Default name for tensorflow layers created. differentiable : boolean - True if the layer should be differentiable, in which case ``odl_op`` - should implement `Operator.derivative` which in turn implements - `Operator.adjoint`, it is properly wrapped in ``tensorflow_layer``, and + ``True`` if the layer should be differentiable, in which case + ``odl_op`` should implement `Operator.derivative` which in turn + implements `Operator.adjoint`. In this case, the adjoint of the + derivative is properly wrapped in ``tensorflow_layer``, and gradients propagate as expected. - Otherwise assumes that the derivative is everywhere zero. + If ``False``, the gradient is defined as everywhere zero. Returns ------- tensorflow_layer : callable - Callable that, when called with an `tensorflow.Tensor` of shape - `(n, *odl_op.domain.shape, 1)` where ``n`` is the number of batches - returns another `tensorflow.Tensor`. + Callable that, when called with a `tensorflow.Tensor` of shape + ``(n, *odl_op.domain.shape, 1)`` where ``n`` is the batch size, + returns another `tensorflow.Tensor` which is a lazy evaluation of + ``odl_op``. If ``odl_op`` is an `Operator`, the shape of the returned tensor is - `(n, *odl_op.range.shape, 1)`. + ``(n, *odl_op.range.shape, 1)``. - Hence for each evaluation, ``odl_op`` is called a total of ``n`` times. + If ``odl_op`` is an `Functional`, the shape of the returned tensor is + ``(n,)``. - The `dtype` of the tensor is the same as the respective ODL spaces. + The ``dtype`` of the tensor is ``odl_op.range.dtype``. """ default_name = name def py_func(func, inp, Tout, stateful=True, name=None, grad=None): """Define custom py_func which takes also a grad op as argument. - We need to overwrite this function since the default tensorflow py_func - does not support custom gradients. + We need to overwrite this function since the default tensorflow + `tf.py_func` does not support custom gradients. + + See tensorflow `issue #1095`_ for more information. Parameters ---------- @@ -71,13 +76,17 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): Input tensors for the function Tout : sequence of `tensorflow.dtype` Datatype of the output(s). - stateful : bool + stateful : bool, optional If the function has internal state, i.e. if calling the function - with a given input always gives the same output. - name : string - Name of the python function - grad : callbable + with a given input repeatedly could give different output. + name : string, optional + Name of the python function. + grad : callbable, optional Gradient of the function. + + References + ---------- + .. _issue #1095: https://github.com/tensorflow/tensorflow/issues/1095 """ if grad is None: return tf.py_func(func, inp, Tout, stateful=stateful, name=name) @@ -97,7 +106,7 @@ def py_func(func, inp, Tout, stateful=True, name=None, grad=None): return tf.py_func(func, inp, Tout, stateful=stateful, name=name) - def tensorflow_layer_grad_impl(x, dx, name): + def tensorflow_layer_grad_impl(x, dy, name): """Implementation of the tensorflow gradient. Gradient in tensorflow is equivalent to the adjoint of the derivative @@ -105,7 +114,7 @@ def tensorflow_layer_grad_impl(x, dx, name): Returns a `tensorflow.Tensor` that represents a lazy application of :: - odl_op.derivative(x).adjoint(dx) + odl_op.derivative(x).adjoint(dy) Parameters ---------- @@ -113,18 +122,18 @@ def tensorflow_layer_grad_impl(x, dx, name): Point(s) in which the derivative should be taken. If ``odl_op`` is an `Operator` the axes are: 0 : batch id. This is a constant if ``fixed_size`` is - true, otherwise it is dynamic. - 1, ..., N-2 : spatial dimensions of data. + ``True``, otherwise it is dynamic. + 1, ..., n-2 : spatial dimensions of data. n-1 : (currently) unused data channel. If ``odl_op`` is a `Functional` the axes are: 0 : batch id. - dx : `tensorflow.Tensor` + dy : `tensorflow.Tensor` Point(s) in which the adjoint of the derivative of the operator should be evaluated. The axes are: 0 : batch id. Should be pairwise matched with ``x``. - 1, ..., M-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. + 1, ..., m-2 : spatial dimensions of data. + m-1 : (currently) unused data channel. name : string Name of the tensor. @@ -134,7 +143,7 @@ def tensorflow_layer_grad_impl(x, dx, name): Lazy result of the computation. If ``odl_op`` is an `Operator` the axes are: 0 : batch id. - 1, ..., N-2 : spatial dimensions of data. + 1, ..., n-2 : spatial dimensions of data. n-1 : (currently) unused data channel. If ``odl_op`` is a `Functional` the axes are: 0 : batch id. @@ -142,7 +151,7 @@ def tensorflow_layer_grad_impl(x, dx, name): with tf.name_scope(name): # Validate the input/output shape x_shape = x.get_shape() - dx_shape = dx.get_shape() + dy_shape = dy.get_shape() try: # Lazy check if the first dimension is dynamic n_x = int(x_shape[0]) @@ -159,16 +168,16 @@ def tensorflow_layer_grad_impl(x, dx, name): assert x_shape[1:] == odl_op.domain.shape + (1,) if odl_op.is_functional: - assert dx_shape[1:] == () + assert dy_shape[1:] == () else: - assert dx_shape[1:] == odl_op.range.shape + (1,) + assert dy_shape[1:] == odl_op.range.shape + (1,) - def _impl(x, dx): + def _impl(x, dy): """Implementation of the adjoint of the derivative. Returns :: - odl_op.derivative(x).adjoint(dx) + odl_op.derivative(x).adjoint(dy) Parameters ---------- @@ -177,17 +186,17 @@ def _impl(x, dx): If ``odl_op`` is an `Operator` the axes are: 0 : batch id. This is a constant if ``fixed_size`` is true, otherwise it is dynamic. - 1, ..., N-2 : spatial dimensions of data. + 1, ..., n-2 : spatial dimensions of data. n-1 : (currently) unused data channel. If ``odl_op`` is a `Functional` the axes are: 0 : batch id. - dx : `numpy.ndarray` + dy : `numpy.ndarray` Point(s) in which the adjoint of the derivative of the operator should be evaluated. The axes are: 0 : batch id. Should be pairwise matched with ``x``. - 1, ..., M-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. + 1, ..., m-2 : spatial dimensions of data. + m-1 : (currently) unused data channel. Returns ------- @@ -196,7 +205,7 @@ def _impl(x, dx): If ``odl_op`` is an `Operator` the axes are: 0 : batch id. - 1, ..., N-2 : spatial dimensions of data. + 1, ..., n-2 : spatial dimensions of data. n-1 : (currently) unused data channel. If ``odl_op`` is a `Functional` the axes are: 0 : batch id. @@ -205,23 +214,23 @@ def _impl(x, dx): if fixed_size: x_out_shape = out_shape assert x.shape == out_shape - assert dx.shape == in_shape + assert dy.shape == in_shape else: x_out_shape = (x.shape[0],) + out_shape[1:] assert x.shape[1:] == out_shape[1:] - assert dx.shape[1:] == in_shape[1:] + assert dy.shape[1:] == in_shape[1:] # Evaluate the operator on all inputs in the batch. out = np.empty(x_out_shape, odl_op.domain.dtype) for i in range(x_out_shape[0]): if odl_op.is_functional: xi = x[i, ..., 0] - dxi = dx[i] - out[i, ..., 0] = np.asarray(odl_op.gradient(xi)) * dxi + dyi = dy[i] + out[i, ..., 0] = np.asarray(odl_op.gradient(xi)) * dyi else: xi = x[i, ..., 0] - dxi = dx[i, ..., 0] - result = odl_op.derivative(xi).adjoint(dxi) + dyi = dy[i, ..., 0] + result = odl_op.derivative(xi).adjoint(dyi) out[i, ..., 0] = np.asarray(result) # Rescale the domain/range according to the weighting since @@ -241,9 +250,9 @@ def _impl(x, dx): return out - with ops.name_scope(name + '_pyfunc', values=[x, dx]) as name_call: + with ops.name_scope(name + '_pyfunc', values=[x, dy]) as name_call: result = py_func(_impl, - [x, dx], + [x, dy], [odl_op.domain.dtype], name=name_call, stateful=False) @@ -266,7 +275,7 @@ def tensorflow_layer(x, name=None): Point(s) to which the layer should be applied. The axes are: 0 : batch id. Can be fixed or dynamic. - 1, ..., M-2 : spatial dimensions of data. + 1, ..., n-2 : spatial dimensions of data. n-1 : (currently) unused data channel. name : string Name of the tensor. Default: Defaultname. @@ -277,8 +286,8 @@ def tensorflow_layer(x, name=None): Lazy result of the computation. If ``odl_op`` is an `Operator` the axes are: 0 : batch id. - 1, ..., N-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. + 1, ..., m-2 : spatial dimensions of data. + m-1 : (currently) unused data channel. If ``odl_op`` is a `Functional` the axes are: 0 : batch id. """ @@ -317,7 +326,7 @@ def _impl(x): The axes are: 0 : batch id. This is a constant if ``fixed_size`` is true, otherwise it is dynamic. - 1, ..., N-2 : spatial dimensions of data. + 1, ..., n-2 : spatial dimensions of data. n-1 : (currently) unused data channel. Returns @@ -326,8 +335,8 @@ def _impl(x): Result of the computation. The axes are: 0 : batch id. Data is pairwise matched with ``x``. - 1, ..., M-2 : spatial dimensions of data. - n-1 : (currently) unused data channel. + 1, ..., m-2 : spatial dimensions of data. + m-1 : (currently) unused data channel. """ # Validate input shape if fixed_size: From a97fcd5ece093f3043c107b277d31afca3d41836 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:48:00 +0200 Subject: [PATCH 68/70] TST: Add rudimentary test for tensorflow --- .../tensorflow/test/tensorflow_test.py | 56 +++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 odl/contrib/tensorflow/test/tensorflow_test.py diff --git a/odl/contrib/tensorflow/test/tensorflow_test.py b/odl/contrib/tensorflow/test/tensorflow_test.py new file mode 100644 index 00000000000..08dff7c0e6a --- /dev/null +++ b/odl/contrib/tensorflow/test/tensorflow_test.py @@ -0,0 +1,56 @@ +# Copyright 2014-2017 The ODL contributors +# +# This file is part of ODL. +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. + +"""Tests for tensorflow.""" + +from __future__ import division +from itertools import permutations +import pytest +import numpy as np +import tensorflow as tf + +import odl +import odl.contrib.tensorflow +from odl.util import all_almost_equal + + +def test_as_tensorflow_layer(): + # Define ODL operator + matrix = np.random.rand(3, 2) + odl_op = odl.MatrixOperator(matrix) + + # Define evaluation points + x = np.random.rand(2) + z = np.random.rand(3) + + # Add empty axes for batch and channel + x_tf = tf.constant(x)[None, ..., None] + z_tf = tf.constant(z)[None, ..., None] + + # Create tensorflow layer from odl operator + odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer( + odl_op, 'MatrixOperator') + y_tf = odl_op_layer(x_tf) + + # Evaluate using tensorflow + result = y_tf.eval().ravel() + expected = odl_op(x) + + assert all_almost_equal(result, expected) + + # Evaluate the adjoint of the derivative, called gradient in tensorflow + result = tf.gradients(y_tf, [x_tf], z_tf)[0].eval().ravel() + expected = odl_op.derivative(x).adjoint(z) + + assert all_almost_equal(result, expected) + + + +if __name__ == '__main__': + with tf.Session(): + pytest.main([str(__file__.replace('\\', '/')), '-v']) From e870da8b5f8a6f405c151c972d4e64392b256686 Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 22 Aug 2017 16:50:32 +0200 Subject: [PATCH 69/70] MAINT: Minor improvements to tensorflow --- odl/contrib/tensorflow/README.md | 2 ++ odl/contrib/tensorflow/space.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/odl/contrib/tensorflow/README.md b/odl/contrib/tensorflow/README.md index f826a4cea7b..dbc02674d11 100644 --- a/odl/contrib/tensorflow/README.md +++ b/odl/contrib/tensorflow/README.md @@ -20,3 +20,5 @@ Specifically * [tensorflow_layer_ray_transform.py](examples/tensorflow_layer_ray_transform.py) shows how a `RayTransform` can be converted to a tensorflow layer. * [tensorflow_operator_matrix.py](examples/tensorflow_operator_matrix.py) shows how `tf.matmul` can be used as a ODL operator. * [tensorflow_tomography.py](examples/tensorflow_tomography.py) shows how tensorflow optimizers can be used with ODL operators to solve inverse problems. + +There are also some rudimentary tests in the [test](test) folder. \ No newline at end of file diff --git a/odl/contrib/tensorflow/space.py b/odl/contrib/tensorflow/space.py index e8811ef79ee..8fc3b719c6a 100644 --- a/odl/contrib/tensorflow/space.py +++ b/odl/contrib/tensorflow/space.py @@ -110,7 +110,7 @@ def data(self, value): self._data = value def __repr__(self): - return '{}.element({})'.format(self.space, self.data) + return '{!r}.element({!r})'.format(self.space, self.data) class TensorflowSpaceOperator(Operator): From 9923fdad526cd503f681be623a650f0f232594da Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 29 Aug 2017 19:18:43 +0200 Subject: [PATCH 70/70] MAINT: Minor style fixes --- odl/contrib/tensorflow/README.md | 6 +++--- odl/tomo/backends/astra_cuda.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/odl/contrib/tensorflow/README.md b/odl/contrib/tensorflow/README.md index dbc02674d11..43ef2864aa0 100644 --- a/odl/contrib/tensorflow/README.md +++ b/odl/contrib/tensorflow/README.md @@ -12,13 +12,13 @@ This package contains ODL functionality related to [TensorFlow](https://www.tens ## Example usage -The [examples](examples) folder contains example on how to use the above functionality. -Specifically +The [examples](examples) folder contains examples on how to use the above functionality. +Specifically: * [tensorflow_layer_matrix.py](examples/tensorflow_layer_matrix.py) shows how an ODL `MatrixOperator` can be converted to a tensorflow layer. * [tensorflow_layer_productspace.py](examples/tensorflow_layer_productspace.py) shows how an ODL operator acting on `ProductSpace`s can be converted to a tensorflow layer. * [tensorflow_layer_ray_transform.py](examples/tensorflow_layer_ray_transform.py) shows how a `RayTransform` can be converted to a tensorflow layer. -* [tensorflow_operator_matrix.py](examples/tensorflow_operator_matrix.py) shows how `tf.matmul` can be used as a ODL operator. +* [tensorflow_operator_matrix.py](examples/tensorflow_operator_matrix.py) shows how `tf.matmul` can be used as an ODL operator. * [tensorflow_tomography.py](examples/tensorflow_tomography.py) shows how tensorflow optimizers can be used with ODL operators to solve inverse problems. There are also some rudimentary tests in the [test](test) folder. \ No newline at end of file diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 63d09d5fab8..2111c17857b 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -65,7 +65,7 @@ def __init__(self, geometry, reco_space, proj_space): # Create a mutually exclusive lock so that two callers cant use the # same shared resource at the same time. - self.mutex = Lock() + self._mutex = Lock() def call_forward(self, vol_data, out=None): """Run an ASTRA forward projection on the given data using the GPU. @@ -84,7 +84,7 @@ def call_forward(self, vol_data, out=None): Projection data resulting from the application of the projector. If ``out`` was provided, the returned object is a reference to it. """ - with self.mutex: + with self._mutex: assert vol_data in self.reco_space if out is not None: assert out in self.proj_space @@ -213,7 +213,7 @@ def __init__(self, geometry, reco_space, proj_space): # Create a mutually exclusive lock so that two callers cant use the # same shared resource at the same time. - self.mutex = Lock() + self._mutex = Lock() def call_backward(self, proj_data, out=None): """Run an ASTRA back-projection on the given data using the GPU. @@ -233,7 +233,7 @@ def call_backward(self, proj_data, out=None): back-projector. If ``out`` was provided, the returned object is a reference to it. """ - with self.mutex: + with self._mutex: assert proj_data in self.proj_space if out is not None: assert out in self.reco_space